#pragma rtGlobals=1		// Use modern global access method.
#include <Multi-peak Fitting 2.0>

Menu "Spectrometry"
	"Iterative Overlayer Correction", iterativeOverLayerGUI()
end

//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//
//																	 //
//																	 //
//  				        	Gold Deconvolution							 //
//																	 //
//																	 //
//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//

function GoldDecon(fwhm,hwhm,boost,iterations,repeats,input_conv_str,energy_Wave_str)

Variable fwhm,hwhm,boost,iterations,repeats
String input_conv_str,energy_Wave_str
Wave input_conv = $input_conv_str, energy_Wave = $energy_Wave_str
Wave input = input_conv

String savDF = getDataFolder(1), new_DataFolder = input_conv_str + "_decon"
NewDataFolder/S/O $new_DataFolder


Variable sigma = (fwhm)/(2*sqrt(2*log(2))), mu = 0, sigma2 = (fwhm)/(2*sqrt(2*log(2)))


Make/O/N=(numpnts(input)) func1, func2, func
Variable scale = ((waveMax(energy_Wave)-waveMin(energy_Wave))/numpnts(input_conv))*numpnts(func1)
KillWaves/Z funcMat
Make/O/N=(numpnts(input),numpnts(input)) funcMat
//Make/O/N=(numpnts(input),numpnts(input)) inputMat
//SetScale x -1,1,funcMat
func1 = exp(-((x-(mu+dimsize(func1,0)/2))^2)/(2*(sigma^2)))/(sigma*sqrt(2*pi))
func2 = (1/pi)*hwhm/(((x+dimsize(func2,0)/2))^2 + hwhm^2)

If(fwhm != 0 && hwhm != 0)
	Duplicate/O func1, func
	Convolve/A func2, func
elseif(fwhm != 0 && hwhm == 0)
	func = func1
elseif(fwhm == 0 && hwhm != 0)
	func = func2
endif
//Resample/UP=(numpnts(func1)*10+9)/DOWN=(numpnts(input)*10) func

funcMat[][] = func[p + numpnts(func)/2 - q]

SetScale x 0,numpnts(func),func, func2
SetScale x 0,numpnts(func),funcMat

Make/O/N=(numpnts(func)) deConv = 1, sumVec

MatrixOp/O toep = funcMat^t x funcMat

Variable p12=0, p13=0, p15=0

Do
	p13 = 0
	Do

			p15 = 0
			Do
				 MatrixOp/O deSum = col(toep, p15) . deConv
			 	sumVec[p15] = deSum[0]
				p15 += 1
			while(p15 < numpnts(deConv))
		
		MatrixOp/O deConv = (input_conv * deConv) / (sumVec)
		p13 += 1
	while(p13 < iterations)
	deConv = (deConv)^(boost)
	p12 += 1
while(p12 < repeats)

Wavestats/Q deConv
deConv = (deConv-V_min)/(V_max-V_min)

MatCon(fwhm,hwhm,deConv,energy_Wave)
KillWaves/Z funcMat, toep
setDataFolder savDF
end

//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//
//																	 //
//																	 //
//  				        	Matrix Convolution							 //
//																	 //
//																	 //
//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//


Function MatCon(fwhm,hwhm,input,energy_Wave)

Variable fwhm,hwhm
Wave input,energy_Wave

Variable sigma = (fwhm)/(2*sqrt(2*log(2))), mu = 0

Make/O/N=(numpnts(input)) func3, func4
Variable scale = ((waveMax(energy_Wave)-waveMin(energy_Wave))/(numpnts(input)))*numpnts(func3)
KillWaves/Z funcMat2
Make/O/N=(numpnts(input),numpnts(input)) funcMat2
Setscale x -scale/2,scale/2,func3, func4
func3 = exp(-((x-mu)^2)/(2*((sigma)^2)))/((sigma)*sqrt(2*pi))
func4 = (1/pi)*(hwhm)/(x^2 + (hwhm)^2)
Duplicate/O func3, funct
Convolve/A func4, funct
//Resample/UP=(numpnts(func3))/DOWN=(numpnts(input)+1) funct

funcMat2[][] = funct[p + numpnts(funct)/2 - q]

MatrixOp/O M_Convolution = funcMat2 x input

WaveStats/Q M_Convolution
M_Convolution = (M_Convolution - V_min)/(V_max-V_min)
KillWaves/Z funcMat2
end

//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//
//																	 //
//																	 //
//  				        		Peak Searching							 //
//																	 //
//																	 //
//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//

//-------------------------------------------------------------------------------------------------------------//
// Second Derivative searching											//
//-------------------------------------------------------------------------------------------------------------//

Function Peak_search(input_wave,energy_wave,wind,repeats)

Wave input_wave,energy_wave
Variable wind, repeats
Make/O/N=(numpnts(input_wave)) first_diff_wave, second_diff_wave, standard_dev

Variable p1 = 0, p2 = 0, p3 = 0 

Do
	first_diff_wave[p1] =  input_wave[min(p1+1,numpnts(input_wave))] - input_wave[min(p1,numpnts(input_wave))]
	second_diff_wave[p1] = input_wave[min(p1+1,numpnts(input_wave))] - 2*input_wave[min(p1,numpnts(input_wave))] + input_wave[max(p1-1,0)]
	standard_dev[p1] =  (input_wave[min(p1+1,numpnts(input_wave))] + 4*input_wave[min(p1,numpnts(input_wave))] + input_wave[max(p1-1,0)])^.5
	p1+=1
while(p1 < numpnts(input_wave))

Do
p1 = 0
	Do
		p2 = 0
		Do
			first_diff_wave[p1] = first_diff_wave[max(min(p1 - wind/2 - 1 +p2,numpnts(input_wave)),0)]
			second_diff_wave[p1] = second_diff_wave[max(min(p1 - wind/2 - 1 +p2,numpnts(input_wave)),0)]
			standard_dev[p1] = standard_dev[max(min(p1 - wind/2 - 1 +p2,numpnts(input_wave)),0)]
			p2+=1
		while(p2 < wind)
		p1+=1
	while(p1 < numpnts(input_wave))
p3+=1
while(p3 < repeats)

Make/O/N=0 peak_position

p1 = 0
Do
	If(abs(first_diff_wave[p1]) < .004 && second_diff_wave[p1] < -0.0048)
		InsertPoints numpnts(peak_position),1,peak_position
		peak_position[numpnts(peak_position) -1] = p1
	endif
p1+=1
while(p1 < numpnts(input_wave))

p2 = 1
//Do
//	if(abs(peak_position[max(p2-1,numpnts(peak_position))] - peak_position[min(p2,numpnts(peak_position))]) <= 2)
//		if(input_wave[max(p2-1,numpnts(peak_position))] < input_wave[min(p2,numpnts(peak_position))])
//			DeletePoints p2-1,1,peak_position
//		else
//			DeletePoints p2,1,peak_position
//		endif
//	endif
//	p2+=1
//while(p2 < numpnts(peak_position))

p3 = 0
Do
	peak_position[p3] = energy_wave[peak_position[p3]]
	p3+=1
while(p3 < numpnts(peak_position))
Return peak_position
end

//-------------------------------------------------------------------------------------------------------------//
// Zero area digital filter peak searching									//
//-------------------------------------------------------------------------------------------------------------//

Function peakFinder(input_wave,energy_wave,chanWidth)

String input_wave,energy_wave
Variable chanWidth
Wave wInput = $input_wave, wEnergy = $energy_wave

Make/O/N=(numpnts(wInput)) filterWave

filterWave[] = 2*wInput[p] - wInput[p-mod(2*chanWidth+1,numpnts(wInput))] - wInput[p+mod(2*chanWidth+1,numpnts(wInput))]

Variable p1=0

Do
	If(filterWave[p1] < .2)
		filterWave[p1] = 0
	endif
	p1+=1
while(p1 < numpnts(wInput))

GoldDecon(.5,.5,1.2,50,20,input_wave,energy_wave)

SetDataFolder filterWave_decon

Wave deConv
Make/O/N=0 peak_position

Variable p2 = 10
Do
	if(deConv[p2] > 0.2)
		InsertPoints numpnts(peak_position),1,peak_position
		peak_position[numpnts(peak_position) -1] = wEnergy[p2]
	endif
	p2+=1
while(p2 < numpnts(deConv))

SetDataFolder ::
Return peak_position
end

//-------------------------------------------------------------------------------------------------------------//
// Zero area digital filter boundry searching								//
//-------------------------------------------------------------------------------------------------------------//

Function zadFilter(input_wave,energy_wave,chanWidth)

String input_wave,energy_wave
Variable chanWidth
Wave wInput = $input_wave, wEnergy = $energy_wave

Make/O/N=(numpnts(wInput)) filterWave

filterWave[] = 2*wInput[p] - wInput[p-mod(2*chanWidth+1,numpnts(wInput))] - wInput[p+mod(2*chanWidth+1,numpnts(wInput))]

Variable p1=0

Do
	If(filterWave[p1] < .2)
		filterWave[p1] = 0
	endif
	p1+=1
while(p1 < numpnts(wInput))

Variable p2 = 0
Make/O/N=(0,2) peakbounds

Do
	If(filterWave[p2-1] < .1 && filterWave[p2] > .1)
		InsertPoints/M=0 numpnts(peakbounds),1,peakbounds
		peakbounds[numpnts(peakbounds) -1][0] = wEnergy[p2]
	endif
	if(filterWave[p2+1] < .1 && filterWave[p2] > .1)
		peakbounds[numpnts(peakbounds) -1][1] = wEnergy[p2]
	endif
	p2+=1
while(p2 < numpnts(wInput))

end

//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//
//																	 //
//																	 //
//  				        	Background Subtraction						 //
//																	 //
//																	 //
//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//

//-------------------------------------------------------------------------------------------------------------//
// LLS operator														//
//-------------------------------------------------------------------------------------------------------------//

Function LLSop(input_wave_str)
String input_wave_str
Wave input_wave = $input_wave_str
String nameHolder = input_wave_str + "_LLS"
Make/O/N=(numpnts(input_wave)) $nameHolder
Wave w1 = $nameHolder

w1[] = ln(ln(sqrt(input_wave[p]+1)+1)+1)

Return w1
end

//-------------------------------------------------------------------------------------------------------------//
// Inverse LLS operator												//
//-------------------------------------------------------------------------------------------------------------//

Function InverseLLSop(input_wave_str)
String input_wave_str
Wave input_wave = $input_wave_str
String nameHolder = input_wave_str + "_invLLS"
Make/O/N=(numpnts(input_wave)) $nameHolder
Wave w1 = $nameHolder

w1[] = ((exp(exp(input_wave[p])-1)-1)^2-1)
Return w1
end

//-------------------------------------------------------------------------------------------------------------//
// Symmetric SNIP procedure											//
//-------------------------------------------------------------------------------------------------------------//

Function SNIP(input_wave_str, order, iterations)
String input_wave_str
Variable order,iterations
Wave input_wave = $input_wave_str
String nameHolder = input_wave_str + "_sub"
Make/O/N=(numpnts(input_wave)) workVect, baseline, LLSop_wave, $nameHolder
Wave w2 = $nameHolder 
Setscale x,pnt2x(input_wave,0),pnt2x(input_wave,numpnts(input_wave)),w2
Make/O/N=5 coefficients
Variable p1 = iterations,p2 = 0
//LLSop_wave = LLSop(input_wave_str)
LLSop_wave[] = input_wave[p]
Wave w1 = LLSop_wave
Do
	p2 = 0
	Do
		coefficients[0] = w1[p2]
		coefficients[1] = (w1[p2-p1]+w1[p2+p1])/2
		coefficients[2] = (-w1[p2-p1]+4*w1[p2-p1/2]+4*w1[p2+p1/2]-w1[p2+p1])/6
		coefficients[3] = (w1[p2-p1]-6*w1[p2-p1*2/3]+15*w1[p2-p1/3]+15*w1[p2+p1/3]-6*w1[p2+p1*2/3]+w1[p2+p1])/20
		coefficients[4] = (-w1[p2-p1]+8*w1[p2-p1*3/4]-28*w1[p2-p1/2]+56*w1[p2-p1/4]+56*w1[p2+p1/4]-28*w1[p2+p1/2]+8*w1[p2+p1*3/4]-w1[p2+p1])/70
		If(order == 2)
			workVect[p2] = min(coefficients[0],coefficients[1])
		elseif(order == 4)
			workVect[p2] = min(coefficients[0],max(coefficients[1],coefficients[2]))
		elseif(order == 6)
			workVect[p2] = min(coefficients[0],max(coefficients[1],max(coefficients[2],coefficients[3])))
		elseif(order == 8)
			workVect[p2] = min(coefficients[0],max(coefficients[1],max(coefficients[2],max(coefficients[3],coefficients[4]))))
		else
			break
		endif
		p2+=1
	while(p2 < numpnts(input_wave))
	w1 = workVect
	p1 -=1
while(p1 > 0)

//baseline = InverseLLSop("LLSop_wave")
baseline[] = w1[p]
WaveStats/Q input_wave
w2 = input_wave - baseline
end

//-------------------------------------------------------------------------------------------------------------//
// Asymmetric SNIP procedure											//
//-------------------------------------------------------------------------------------------------------------//

Function aSNIP(input_wave_str, order, iterationsR,iterationsL)
String input_wave_str
Variable order,iterationsR,iterationsL
Wave input_wave = $input_wave_str
String nameHolder = input_wave_str + "_sub"
KillWaves/Z workVect, baseline, LLSop_wave, $nameHolder,coefficients
Duplicate/O input_wave, workVect, baseline, LLSop_wave, $nameHolder
Wave w2 = $nameHolder 
Make/O/N=5 coefficients
Variable p1 = iterationsR,p2 = 0, p3 = iterationsL
//LLSop_wave = LLSop(input_wave_str)
LLSop_wave[] = input_wave[p]
Wave w1 = LLSop_wave
Do
	p2 = 0
	Do
		coefficients[0] = w1[p2]
		coefficients[1] = (w1[p2-p1]+w1[p2+p3])/2
		coefficients[2] = (-w1[p2-p1]+4*w1[p2-p1/2]+4*w1[p2+p3/2]-w1[p2+p3])/6
		coefficients[3] = (w1[p2-p1]-6*w1[p2-p1*2/3]+15*w1[p2-p1/3]+15*w1[p2+p3/3]-6*w1[p2+p3*2/3]+w1[p2+p3])/20
		coefficients[4] = (-w1[p2-p1]+8*w1[p2-p1*3/4]-28*w1[p2-p1/2]+56*w1[p2-p1/4]+56*w1[p2+p3/4]-28*w1[p2+p3/2]+8*w1[p2+p3*3/4]-w1[p2+p3])/70
		If(order == 2)
			workVect[p2] = min(coefficients[0],coefficients[1])
		elseif(order == 4)
			workVect[p2] = min(coefficients[0],max(coefficients[1],coefficients[2]))
		elseif(order == 6)
			workVect[p2] = min(coefficients[0],max(coefficients[1],max(coefficients[2],coefficients[3])))
		elseif(order == 8)
			workVect[p2] = min(coefficients[0],max(coefficients[1],max(coefficients[2],max(coefficients[3],coefficients[4]))))
		else
			break
		endif
		p2+=1
	while(p2 <= numpnts(input_wave))
	w1[] = workVect[p]
	if(p1 > 2)
		p1 -=1
	endif
	if(p3 > 2)
		p3 -=1
	endif
while(p1 > 2 && p3 > 2)

//baseline = InverseLLSop("LLSop_wave")
baseline[] = w1[p]
w2 = input_wave - baseline
end


//-------------------------------------------------------------------------------------------------------------//
// Scaling procedure													//
//-------------------------------------------------------------------------------------------------------------//

Function scaleWaves(wList)
String wList
String nameHolder
Make/O/N=(itemsinlist(wList)) wStats = 0
Variable p1 = 0, p2 = 0

Do
	nameholder = stringfromlist(p1,wList)
	Wavestats/Q/R=[0,30] $nameholder
	wStats[p1] = V_adev 
	p1+=1
while(p1 < itemsinlist(wList))

Variable minwStats = Wavemin(wStats)
wStats = minwStats/wStats
Do
	nameHolder = stringfromlist(p2,wList) + "_scaled"
	Make/O/N=(numpnts($stringfromlist(p2,wList))) $nameHolder
	Wave w2  =$stringfromlist(p2,wList), w1 = $nameHolder
	w1 = w2*wStats[p2]
	p2+=1
while(p2 < numpnts(wStats))

end


//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//
//																	 //
//																	 //
//  				        		Batch Peak Fitting						 //
//																	 //
//																	 //
//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//

Function batchFit(ywList,xwList,type)
String ywList,xwList,type

MPF2_AutoMPFit("PeakResults", type, "PeakCoefs_%d", "None", "Background", ywList, xwList, 5)

end


//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//
//																	 //
//																	 //
//  				        Roughness Penalty Approach						 //
//																	 //
//																	 //
//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//

Function rPSmth(inputWave, [alphaMax, alphaMin, stepSize])
Wave inputWave
Variable alphaMax,alphaMin,stepSize

String wName = NameOfWave(inputWave)
Duplicate/O inputWave,$(wName + "_smth")
Wave w1 =  $(wName + "_smth")

If(alphaMax == 0)
	alphaMax = 10
endif

If(stepSize == 0)
	stepSize = .1
endif

Make/O/N=(dimsize(inputWave,0),(dimsize(inputWave,0)-2)) qMatrix
Make/O/N=(dimsize(inputWave,0)-2,dimsize(inputWave,0)-2) rMatrix

Variable p1 = 0, p2 = 0

Do
	qMatrix[p2-1][p2] = 1/(p2+1-p2)
	qMatrix[p2][p2] = -1/(p2+1-p2) - 1/(p2+2 - p2+1)
	qMatrix[p2+1][p2] = 1/(p2+2 - p2+1)
	p2 += 1
while(p2 < dimsize(qMatrix, 1))

p2 = 0
Do
	rMatrix[p2][p2+1] =  1/6*(p2+2 - p2+1)
	rMatrix[p2][p2] = 1/3*((p2+1-p2) + (p2+2 - p2+1))
	rMatrix[p2+1][p2] = 1/6*(p2+2 - p2+1)
	p2 += 1
while(p2 < dimsize(rMatrix, 1))

MatrixOp/O kMatrix = qMatrix x inv(rMatrix) x qMatrix^t

Make/O/N=(numpnts(inputWave),numpnts(inputWave)) identityMatrix

p2 = 0
Do
	identityMatrix[p2][p2] = 1
	p2 += 1
while(p2 < dimsize(identityMatrix, 1))
//MatrixOp/O aMatrix = inv(identityMatrix + Alpha*kMatrix)
//MatrixOp/O smoothed = aMatrix x inputWave



Variable p3 = 0, alpha, alphaPnts = (alphaMax-alphaMin)/stepSize

if(alphaMax != alphaMin)
	KillWaves/Z GCV
	Make/O/N=(alphaPnts) GCV

	If(alphaMin != 0)
		Setscale x,alphaMin,alphaMax,GCV
	else
		Setscale x,alphaMin + stepSize,alphaMax + stepsize,GCV
	endif

	Do
		alpha = pnt2x(GCV,p3)
		MatrixOp/O aMatrix = inv(identityMatrix + alpha*kMatrix)
		MatrixOp/O difference = (inputWave - aMatrix x inputWave)
		difference = difference^2
		GCV[p3] = sum(difference)/(numpnts(inputWave)*(1-1/numpnts(inputWave)*(MatrixTrace(aMatrix)))^2)
		p3 += 1
	while(p3 < numpnts(GCV))
	WaveStats/Q GCV
	alpha = pnt2x(GCV, V_minRowLoc)
else
	alpha = alphaMax
endif
MatrixOp/O aMatrix = inv(identityMatrix + alpha*kMatrix)
MatrixOp/O w1= aMatrix x inputWave
MatrixOp/O edfMatrix = identityMatrix - aMatrix
Variable/G EDF = MatrixTrace(edfMatrix)
Killwaves/Z aMatrix, edfMatrix, identityMatrix, kMatrix, difference, qMatrix, rMatrix
//MatrixOp/O secondDeriv_Smooth = inv(rMatrix + alpha*(qMatrix^t x qMatrix)) x qMatrix^t x inputWave
end

//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//
//																	 //
//																	 //
//  			            UPS Secondary Electron Subtraction				 //
//																	 //
//																	 //
//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//

Function secondaryElecSub(inputSpec)
Wave inputSpec
String inputSpecName = NameOfWave(inputSpec)

end

//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//
//																	 //
//																	 //
//  					Iterative Carbon Overlayer Correction				 //
//																	 //
//																	 //
//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//

Function iterativeOverlayerCorrection(peakAreas, atomicMasses, elementList, carbonArea, theta, [distance, density, AtomicWeight, BandGap, numValanceElectrons, photonEnergy])
Wave peakAreas, atomicMasses
String elementList
Variable carbonArea, theta, distance, density, AtomicWeight, BandGap, numValanceElectrons, photonEnergy
Variable mValue = .474341649025, densityCarbon = 2.25*.6022, atomicMassCarbon = 12
theta = pi/180*theta
if(ParamIsDefault(distance))
	distance = 15
endif
if(ParamIsDefault(density))
	density = 6.2
endif
if(ParamIsDefault(AtomicWeight))
	AtomicWeight = 546.458
endif
if(ParamIsDefault(BandGap))
	BandGap = 3.05
endif
if(ParamIsDefault(numValanceElectrons))
	numValanceElectrons = 66
endif
if(ParamIsDefault(photonEnergy))
	photonEnergy = 1486.7
endif

if(!WaveExists($"root:crossSections:C1s"))
	getCrossSections(root:SensitivityFactors:C1s)
endif
Variable kPrime, kValue, carbonAreaPredicted = 0, distanceTemp = distance, distanceStep = 1, direction = 1
Print "Density: ", density
Print "AtomicWeight:", AtomicWeight
Print "BandGap: ", BandGap
Print "numValanceElectrons: ", numValanceElectrons
Make/D/O/N=(numpnts(peakAreas)) kPrimeNValues, nValues, nAiValues
density *= .6022
Do
	Variable p1,p2
	for(p1 = 0; p1 < numpnts(kPrimeNValues); p1 +=1)
		Wave waveTemp = $("root:SensitivityFactors:" + stringFromList(p1, elementList))
		if(!WaveExists($("root:crossSections:" + stringFromList(p1, elementList))))
			getCrossSections(waveTemp)
		endif
		Wave waveTemp1 = $("root:crossSections:" + stringFromList(p1, elementList))           
		Variable IMFP = TPP2M(waveTemp[0], density = density/.6022, AtomicWeight = AtomicWeight, BandGap = BandGap, numValanceElectrons =  numValanceElectrons, photonEnergy = photonEnergy)
		//Print cos(theta)
		//Print "IMFP: ", IMFP
		//Print "IMFP Carbon: ", mValue*(photonEnergy - waveTemp[0])^(1/2)
		//(1+(1/2)*waveTemp1[4]*((3/2)*sin((pi/180)*45)^2-1))
		kPrimeNValues[p1] = peakAreas[p1]*spinSplitAreaFactor(stringFromList(p1,elementList))/(waveTemp1[1]*.01*IMFP*exp(- distanceTemp/(mValue*(photonEnergy - waveTemp[0])^(1/2)*cos(pi/2 - theta))))
	endfor
	kPrime = sum(kPrimeNValues)
	nValues = kPrimeNValues/kPrime
	nAiValues = nValues*atomicMasses
	
	//for(p1 = 0; p1 < numpnts(kPrimeNValues); p1 +=1)
	//	Wave waveTemp = $("root:SensitivityFactors:" + stringFromList(p1, elementList))
	//	IMFP = TPP2M(waveTemp[0], density = density, AtomicWeight = AtomicWeight, BandGap = BandGap, numValanceElectrons = numValanceElectrons, photonEnergy = photonEnergy)
		//Print ("Area " + stringFromList(p1, elementList) + ":"), kValue*nValues[p1]*getCrossSections(waveTemp)*IMFP*exp(- distanceTemp/(mValue*(photonEnergy - waveTemp[0])^(1/2)*cos(theta)))
		//Print "IMFP: ", IMFP
//	endfor
	
	Wave wCarbon = root:crossSections:C1s
	kValue = kPrime*sum(nAiValues)/density
	//(1+(1/2)*waveTemp1[4]*((3/2)*sin((pi/180)*45)^2-1))
	carbonAreaPredicted = kValue*densityCarbon/(atomicMassCarbon)*wCarbon[1]*.01*mValue*(photonEnergy - 284.2)^(1/2)*(1 - exp(- distanceTemp/(mValue*(photonEnergy - 284.2)^(1/2)*cos(pi/2 - theta))))
	if(carbonArea - carbonAreaPredicted < 0)
		if(direction > 0)
			direction = -abs(direction)
			distanceStep /= 10
		endif
		distanceTemp -= distanceStep
	else
		if(direction < 0)
			direction = abs(direction)
			distanceStep /= 10
		endif
		distanceTemp += distanceStep
	endif
	p2 += 1
While(abs(carbonArea - carbonAreaPredicted) > .01)
Print carbonArea, carbonAreaPredicted, p2
Print "Distance: ", distanceTemp
end

Function TPP2M(energy, [density,AtomicWeight,BandGap,numValanceElectrons, photonEnergy])
Variable energy,density,AtomicWeight,BandGap,numValanceElectrons, photonEnergy
if(ParamIsDefault(density))
	density = 6.2
endif
if(ParamIsDefault(AtomicWeight))
	AtomicWeight = 546.458
endif
if(ParamIsDefault(BandGap))
	BandGap = 3.05
endif
if(ParamIsDefault(numValanceElectrons))
	numValanceElectrons = 66
endif
if(ParamIsDefault(photonEnergy))
	photonEnergy = 1486.7
endif
Variable Ep = 28.8*(density*numValanceElectrons/AtomicWeight)^(1/2)
Variable betaVar = -0.1 + 0.944*(Ep^2 + BandGap^2)^(-1/2) + .069*density^(0.1)
Variable gammaVar = 0.191*density^(-1/2)
Variable uVar = density*numValanceElectrons/AtomicWeight
Variable cVar = 1.97 -0.91*uVar
Variable dVar = 53.4 - 20.8*uVar
Return (photonEnergy - energy)/(Ep^2 * (betaVar*ln(gammaVar*(photonEnergy - energy)) - cVar/(photonEnergy - energy) + dVar/(photonEnergy - energy)^2))
end

Function attenuationLength(energy, [zNumber,density,AtomicWeight,BandGap,numValanceElectrons, photonEnergy])
Variable energy,zNumber,density,AtomicWeight,BandGap,numValanceElectrons,photonEnergy
if(ParamIsDefault(zNumber))
	zNumber = 1
endif
if(ParamIsDefault(density))
	density = 6.2
endif
if(ParamIsDefault(AtomicWeight))
	AtomicWeight = 546.45
endif
if(ParamIsDefault(BandGap))
	BandGap = 3.05
endif
if(ParamIsDefault(numValanceElectrons))
	numValanceElectrons = 160
endif
if(ParamIsDefault(photonEnergy))
	photonEnergy = 1486.7
endif
	Return (1 - 0.028*zNumber^(1/2))*(0.501 + 0.068*ln(photonEnergy - energy))*TPP2M(energy, density = density, AtomicWeight = AtomicWeight, BandGap = BandGap, numValanceElectrons = numValanceElectrons, photonEnergy = photonEnergy)
end

Function getCrossSections(element)
Wave element
String savDF = GetDataFolder(1)
NewDataFolder/O/S root:crossSections
String fixedName
Variable tempNum = str2num(NameOfWave(element)[strLen(NameOfWave(element))-1,inf])
if(tempNum)
	fixedName = LowerStr(NameOfWave(element))[0,strlen(NameOfWave(element))-2]
else
	fixedName = LowerStr(NameOfWave(element))[0,strlen(NameOfWave(element))]
endif
String urlString = "http://ulisse.elettra.trieste.it/services/elements/data/" + fixedName + ".txt"
Variable photonEnergy, crossSection1, crossSection2, crossSection3, asym1, asym2, asym3
String crossSections = FetchURL(urlString)
sscanf crossSections[strsearch(crossSections, "1486", 0), strsearch(crossSections, "\n", strsearch(crossSections, "1486", 0))], "%f %f %f %f %f %f %f", photonEnergy, crossSection1, crossSection2, crossSection3, asym1, asym2, asym3
if(photonEnergy)
Make/O/N=7 $NameOfWave(element) = {photonEnergy, crossSection1, crossSection2, crossSection3, asym1, asym2, asym3}
endif
SetDataFolder savDF
if(photonEnergy)
Wave w1 = $("root:crossSections:" + NameOfWave(element))
Return w1[2]
else
Return 0
endif
end

Function defineSensitivityFactors()
String savDF = GetDataFolder(1)
NewDataFolder/O/S root:SensitivityFactors
Make/O Al2p = {74,0.256,0.249,88,68,87,69,3,1,2}
Make/O Al2s = {119,0.312,0.230,133,113,132,114,2,1,2}
Make/O Sb3d3 = {539,3.649,3.318,549,529,548,530,2,1,2}
Make/O Sb3d5 = {530,5.334,5.334,543,523,542,524,3,1,2}
Make/O Sb3p1 = {814,0.809,0.355,824,804,823,805,2,1,2}
Make/O Sb3p3 = {768,1.733,0.657,778,758,777,759,2,1,2}
Make/O Sb4d = {34,1.405,1.208,48,28,47,29,2,1,2}
Make/O Sb_MNN = {1022,0.000,0.883,1040,1010,1039,1011,2,1,2}
Make/O Sb_MNN1 = {1032,0.000,0.866,1040,1010,1039,1011,0,1,2}
Make/O Ar2p = {241,1.231,0.962,255,235,254,236,3,1,2}
Make/O Ar2s = {319,0.494,0.321,332,312,331,313,0,1,2}
Make/O Ar_LMM = {1272,0.000,0.266,1285,1260,1284,1261,0,1,2}
Make/O As2p1 = {1362,1.232,0.069,1372,1352,1371,1353,2,1,2}
Make/O As2p3 = {1326,2.464,0.525,1340,1320,1339,1321,2,1,2}
Make/O As3d = {44,0.741,0.713,56,36,55,37,3,1,2}
Make/O As3p = {142,1.341,0.577,162,137,161,138,0,1,2}
Make/O As3s = {205,0.350,0.577,215,195,214,196,0,1,2}
Make/O As_LMM = {266,0.000,0.340,284,259,283,260,2,1,2}
Make/O As_LMM1 = {230,0.000,0.628,250,220,249,221,0,1,2}
Make/O Ba3d3 = {795,5.001,4.835,805,785,804,786,2,1,2}
Make/O Ba3d5 = {780,7.343,6.971,794,774,793,775,3,1,2}
Make/O Ba4d = {90,3.224,1.852,105,85,104,86,2,1,2}
Make/O Ba4p3 = {177,0.620,0.554,187,167,186,168,0,1,2}
Make/O Ba_MNN = {886,0.000,0.682,930,860,929,861,2,1,2}
Make/O Ba_MNN1 = {900,0.000,0.962,930,860,929,861,0,1,2}
Make/O Be1s = {113,0.080,0.080,127,107,126,108,3,1,2}
Make/O Bi4d3 = {467,1.952,0.792,477,457,476,458,2,1,2}
Make/O Bi4d5 = {443,2.907,1.163,462,437,461,438,2,1,2}
Make/O Bi4f = {159,9.848,5.656,172,152,171,153,3,1,2}
Make/O Bi4f5 = {164,4.515,4.515,174,154,173,155,2,1,2}
Make/O Bi4f7 = {159,5.656,5.656,178,153,177,154,2,1,2}
Make/O Bi4p3 = {681,1.606,0.412,691,671,690,672,2,1,2}
Make/O Bi5d = {26,1.551,0.916,40,20,39,21,0,1,2}
Make/O B1s = {191,0.171,0.171,202,182,201,183,3,1,2}
Make/O B_KLL = {1310,0.000,0.006,1325,1290,1324,1291,0,1,2}
Make/O Br3d = {69,1.149,1.066,83,63,82,64,3,1,2}
Make/O Br3p = {182,1.502,0.711,198,178,197,179,2,1,2}
Make/O Br3p1 = {189,0.512,0.289,199,179,198,180,0,1,2}
Make/O Br3p3 = {182,0.988,0.711,192,172,191,173,0,1,2}
Make/O Br3s = {256,0.410,0.178,266,246,265,247,0,1,2}
Make/O Cd3d3 = {412,2.858,2.858,422,402,421,403,2,1,2}
Make/O Cd3d5 = {405,4.148,4.148,419,399,418,400,3,1,2}
Make/O Cd3p1 = {652,0.820,0.338,662,642,661,643,2,1,2}
Make/O Cd3p3 = {618,1.682,0.788,633,613,632,614,2,1,2}
Make/O Cd3s = {772,0.800,0.788,782,762,781,763,0,1,2}
Make/O Cd4d = {11,0.752,0.582,21,1,20,2,2,1,2}
Make/O Cd_MNN = {1103,0.000,0.673,1120,1090,1119,1091,2,1,2}
Make/O Cd_MNN1 = {1110,0.000,0.938,1120,1090,1119,1091,0,1,2}
Make/O Ca2p = {347,1.927,1.281,361,341,360,342,3,1,2}
Make/O Ca2p1 = {350,0.658,0.585,360,340,359,341,2,1,2}
Make/O Ca2p3 = {347,1.281,1.281,362,342,361,343,2,1,2}
Make/O Ca2s = {439,0.548,0.303,454,434,453,435,2,1,2}
Make/O Ca_LMM = {1197,0.000,0.212,1207,1175,1206,1176,0,1,2}
Make/O C1s = {284.8,0.314,0.314,298,278,297,279,3,1,2}
Make/O C_KLL = {1223,0.000,0.020,1240,1200,1239,1201,2,1,2}
Make/O Ce3d = {882,12.000,1.560,916,876,915,877,2,1,2}
Make/O Ce3d3 = {900,5.001,1.052,910,890,909,891,2,1,2}
Make/O Ce3d5 = {882,7.343,1.560,902,872,901,873,3,1,2}
Make/O Ce4d = {108,2.720,0.817,125,105,124,106,2,1,2}
Make/O Ce_MNN = {833,0.000,0.838,860,820,859,821,2,1,2}
Make/O Cs3d = {724,11.800,6.801,759,719,758,720,3,1,2}
Make/O Cs3d3 = {738,4.805,4.709,748,728,747,729,2,1,2}
Make/O Cs3d5 = {724,6.995,6.801,738,718,737,719,2,1,2}
Make/O Cs3p3 = {997,1.658,0.460,1007,987,1006,988,2,1,2}
Make/O Cs4d = {75,2.761,1.726,90,70,89,71,2,1,2}
Make/O Cs4p3 = {158,0.620,0.533,168,148,167,149,0,1,2}
Make/O Cs_MNN = {918,0.000,1.092,950,890,949,891,2,1,2}
Make/O Cs_MNN1 = {931,0.000,0.976,950,890,949,891,0,1,2}
Make/O Cl2p = {199,0.954,0.797,213,193,212,194,3,1,2}
Make/O Cl2s = {270,0.468,0.316,285,265,284,266,2,1,2}
Make/O Cl_LMM = {1304,0.000,0.181,1315,1290,1314,1291,0,1,2}
Make/O Cr2p = {577,2.488,1.623,607,567,606,568,0,1,2}
Make/O Cr2p1 = {586,0.860,0.860,596,576,595,576,0,1,2}
Make/O Cr2p3 = {577,1.623,1.623,589,569,588,570,3,1,2}
Make/O Cr2s = {698,0.892,0.109,708,688,707,689,0,1,2}
Make/O Cr3p = {45,0.200,0.238,58,38,57,39,2,1,2}
Make/O Cr3s = {77,0.080,0.069,87,67,86,68,0,1,2}
Make/O Cr_LMM = {959,0.000,0.118,975,945,974,946,0,1,2}
Make/O Cr_LMM1 = {997,0.000,0.361,1015,980,1014,981,0,1,2}
Make/O Cr_LMM2 = {917,0.000,0.366,925,905,924,906,0,1,2}
Make/O Co2p = {781,3.529,2.322,815,765,814,766,0,1,2}
Make/O Co2p1 = {796,1.056,1.008,806,786,805,787,0,1,2}
Make/O Co2p3 = {781,2.113,2.322,795,770,794,771,3,1,2}
Make/O Co3p = {61,0.486,0.347,75,55,74,56,0,1,2}
Make/O Co3s = {103,0.170,0.117,113,93,112,94,0,1,2}
Make/O Co_LMM = {713,0.000,0.000,735,685,734,686,2,1,2}
Make/O Co_LMM1 = {777,0.000,0.438,795,765,794,766,0,1,2}
Make/O Co_LMM2 = {838,0.000,0.302,865,820,864,821,0,1,2}
Make/O Cu2p = {934,4.395,2.930,974,924,973,925,0,1,2}
Make/O Cu2p1 = {954,1.313,1.152,964,944,963,945,0,1,2}
Make/O Cu2p3 = {934,2.626,2.930,952,927,951,928,3,1,2}
Make/O Cu3p = {77,0.897,0.552,109,59,109,60,2,1,2}
Make/O Cu_LMM = {568,0.000,0.129,620,540,619,541,2,1,2}
Make/O Cu_LMM1 = {548,0.000,1.087,560,530,559,531,0,1,2}
Make/O Cu_LMM2 = {640,0.000,0.134,660,620,659,621,0,1,2}
Make/O Cu_LMM3 = {648,0.000,0.189,660,620,659,621,0,1,2}
Make/O Dy3d5 = {1295,7.343,0.000,1305,1285,1304,1286,2,1,2}
Make/O Dy4d = {154,2.668,2.668,175,145,174,146,3,1,2}
Make/O Dy4p3 = {293,0.750,0.751,306,286,305,287,2,1,2}
Make/O Er4d = {169,2.650,2.650,180,160,179,161,3,1,2}
Make/O Er4p = {320,0.741,0.741,335,315,334,316,2,1,2}
Make/O Eu3d = {1126,12.000,1.890,1185,1110,1184,1111,2,1,2}
Make/O Eu3d5 = {1126,7.343,1.890,1150,1115,1149,1116,3,1,2}
Make/O Eu4d = {128,2.689,2.689,180,120,179,121,2,1,2}
Make/O F1s = {686,1.000,1.000,699,679,698,680,3,1,2}
Make/O F2s = {30,0.056,0.056,42,22,41,23,2,1,2}
Make/O F_KLL = {832,0.000,0.430,845,820,844,821,2,1,2}
Make/O F_KLL1 = {858,0.000,0.029,870,845,869,846,0,1,2}
Make/O F_KLL2 = {877,0.000,0.028,890,870,889,871,0,1,2}
Make/O Gd3d = {1186,12.000,0.896,1215,1175,1214,1176,2,1,2}
Make/O Gd4d = {141,2.683,2.683,167,132,166,133,3,1,2}
Make/O Ga2p1 = {1144,1.376,1.023,1154,1134,1153,1135,2,1,2}
Make/O Ga2p3 = {1117,2.751,2.203,1140,1110,1139,1111,3,1,2}
Make/O Ga3d = {20,0.438,0.438,34,14,33,15,2,1,2}
Make/O Ga3p = {105,1.145,0.545,121,101,120,102,0,1,2}
Make/O Ga3s = {160,0.290,0.438,170,150,169,151,0,1,2}
Make/O Ga_LMM = {419,0.000,1.008,430,410,429,411,0,1,2}
Make/O Ga_LMM1 = {392,0.000,1.991,405,380,404,381,0,1,2}
Make/O Ga_LMM2 = {487,0.000,0.241,500,460,499,461,0,1,2}
Make/O Ga_LMM3 = {504,0.000,0.239,515,490,514,491,0,1,2}
Make/O Ga_LMM4 = {514,0.000,0.356,520,505,519,506,0,1,2}
Make/O Ge2p1 = {1250,1.368,0.711,1260,1240,1259,1241,2,1,2}
Make/O Ge2p3 = {1219,2.735,1.479,1245,1210,1244,1211,3,1,2}
Make/O Ge3d = {31,0.535,0.520,45,25,44,26,2,1,2}
Make/O Ge3p = {124,1.243,0.541,140,120,139,121,0,1,2}
Make/O Ge3s = {181,0.310,0.541,191,171,190,172,0,1,2}
Make/O Ge_LMM = {342,0.000,0.915,350,330,349,331,2,1,2}
Make/O Ge_LMM1 = {310,0.000,1.934,320,300,319,301,0,1,2}
Make/O Ge_LMM2 = {412,0.000,0.249,422,402,421,403,0,1,2}
Make/O Ge_LMM3 = {433,0.000,0.247,443,413,442,414,0,1,2}
Make/O Ge_LMM4 = {444,0.000,0.369,464,424,463,425,0,1,2}
Make/O Au4d3 = {354,1.702,0.790,364,344,363,345,2,1,2}
Make/O Au4d5 = {336,2.513,1.103,355,335,354,336,2,1,2}
Make/O Au4f = {85,6.805,3.849,99,79,98,80,3,1,2}
Make/O Au4f5 = {89,3.019,3.019,99,79,98,80,2,1,2}
Make/O Au4f7 = {85,3.849,3.849,100,80,99,81,2,1,2}
Make/O Au4p1 = {643,0.382,0.072,653,633,652,634,2,1,2}
Make/O Au4p3 = {547,1.134,0.363,557,537,556,538,2,1,2}
Make/O Hf4d = {213,3.052,1.169,243,203,242,204,0,1,2}
Make/O Hf4d3 = {224,1.292,0.853,234,214,233,215,2,1,2}
Make/O Hf4d5 = {213,1.844,1.169,243,203,242,204,2,1,2}
Make/O Hf4f = {17,2.901,2.405,29,9,28,10,3,1,2}
Make/O Hf4p1 = {437,0.385,0.187,447,427,446,428,0,1,2}
Make/O Hf4p3 = {382,0.810,0.456,392,372,391,373,2,1,2}
Make/O Ho4d = {161,2.659,2.659,185,155,184,156,3,1,2}
Make/O Ho4p = {306,0.746,0.746,320,300,319,301,2,1,2}
Make/O In3d3 = {453,3.239,3.239,463,443,462,444,2,1,2}
Make/O In3d5 = {445,4.530,4.530,459,439,458,440,3,1,2}
Make/O In3p1 = {704,0.829,0.365,714,694,713,695,2,1,2}
Make/O In3p3 = {666,1.705,0.761,680,660,679,661,2,1,2}
Make/O In3s = {828,0.830,0.761,838,818,837,819,0,1,2}
Make/O In4d = {17,0.840,0.761,27,7,26,8,2,1,2}
Make/O In_MNN = {1076,0.000,0.784,1090,1060,1089,1061,2,1,2}
Make/O In_MNN1 = {1084,0.000,0.769,1090,1060,1089,1061,0,1,2}
Make/O I3d3 = {630,4.168,4.168,640,620,639,621,2,1,2}
Make/O I3d5 = {619,6.302,6.302,633,613,632,614,3,1,2}
Make/O I3p1 = {930,0.779,0.211,940,920,939,921,2,1,2}
Make/O I3p3 = {874,1.700,0.476,884,864,883,865,2,1,2}
Make/O I4d = {50,2.009,1.507,64,44,63,45,2,1,2}
Make/O I_MNN = {982,0.000,0.784,1000,950,999,951,2,1,2}
Make/O I_MNN1 = {971,0.000,0.808,1000,950,999,951,0,1,2}
Make/O Ir4d3 = {313,1.518,0.744,323,303,322,304,2,1,2}
Make/O Ir4d5 = {297,2.298,1.124,313,288,312,289,2,1,2}
Make/O Ir4f = {62,5.484,3.332,75,55,74,56,3,1,2}
Make/O Ir4f5 = {65,2.286,2.496,75,55,74,56,2,1,2}
Make/O Ir4f7 = {62,3.124,3.332,74,54,73,55,2,1,2}
Make/O Ir4p1 = {579,0.389,0.151,589,569,588,570,2,1,2}
Make/O Ir4p3 = {497,1.063,0.577,507,487,506,488,2,1,2}
Make/O Fe2p = {710,2.946,1.964,740,700,739,701,0,1,2}
Make/O Fe2p1 = {723,0.972,0.972,733,713,732,714,0,1,2}
Make/O Fe2p3 = {710,1.964,1.964,727,702,726,703,3,1,2}
Make/O Fe3p = {55,0.362,0.292,68,48,67,49,2,1,2}
Make/O Fe3s = {92,0.140,0.243,102,82,101,83,0,1,2}
Make/O Fe_LMM = {784,0.000,0.506,795,770,794,771,2,1,2}
Make/O Fe_LMM1 = {839,0.000,0.486,850,820,849,821,0,1,2}
Make/O Fe_LMM2 = {888,0.000,0.401,900,850,899,851,0,1,2}
Make/O Kr3d = {88,1.400,1.249,102,82,101,83,3,1,2}
Make/O Kr3p = {208,1.601,0.781,230,205,229,206,0,1,2}
Make/O Kr3p1 = {216,0.533,0.781,230,205,229,206,0,1,2}
Make/O Kr3p3 = {208,1.068,0.781,230,205,229,206,0,1,2}
Make/O La3d = {834,12.000,5.135,867,827,866,828,2,1,2}
Make/O La3d3 = {851,5.001,3.316,861,841,860,842,2,1,2}
Make/O La3d5 = {834,7.343,5.135,850,825,849,826,3,1,2}
Make/O La4d = {101,2.730,2.730,115,95,114,96,2,1,2}
Make/O La4p3 = {195,0.650,0.890,205,185,204,186,0,1,2}
Make/O La_MNN = {867,0.000,1.907,880,840,879,841,2,1,2}
Make/O Pb4d3 = {435,1.868,0.957,445,425,444,426,2,1,2}
Make/O Pb4d5 = {413,2.774,1.180,432,407,431,408,2,1,2}
Make/O Pb4f = {138,9.000,5.131,152,132,151,133,3,1,2}
Make/O Pb4f5 = {143,3.820,4.289,153,133,152,134,2,1,2}
Make/O Pb4f7 = {138,5.172,5.131,150,125,149,126,2,1,2}
Make/O Pb4p3 = {645,1.175,0.392,655,635,654,636,2,1,2}
Make/O Pb5d = {19,1.414,0.848,34,14,33,15,0,1,2}
Make/O Li1s = {56,0.028,0.028,71,51,70,52,3,1,2}
Make/O Lu4d = {197,2.617,2.617,218,193,217,194,2,1,2}
Make/O Lu4f = {7,0.000,0.000,20,0,19,1,3,1,2}
Make/O Lu4p1 = {413,0.385,0.727,423,403,422,404,2,1,2}
Make/O Lu4p3 = {359,0.800,0.727,374,354,373,355,2,1,2}
Make/O Lu4s = {509,0.420,0.727,519,499,518,500,0,1,2}
Make/O Mg1s = {1305,1.035,0.976,1317,1297,1316,1298,3,1,2}
Make/O Mg2p = {51,0.167,0.167,64,44,63,45,2,1,2}
Make/O Mg2s = {90,0.274,0.206,104,84,103,85,2,1,2}
Make/O Mg_KLL = {305,0.000,2.627,318,298,317,299,2,1,2}
Make/O Mn2p = {641,2.688,1.757,672,632,671,633,0,1,2}
Make/O Mn2p1 = {652,0.923,0.923,662,642,661,643,0,1,2}
Make/O Mn2p3 = {641,1.757,1.757,654,634,653,635,3,1,2}
Make/O Mn3p = {48,0.260,0.265,63,43,62,44,2,1,2}
Make/O Mn3s = {83,0.110,0.219,93,53,92,54,0,1,2}
Make/O Mn_LMM = {900,0.000,0.347,920,890,919,891,2,1,2}
Make/O Mn_LMM1 = {852,0.000,0.648,870,840,869,841,0,1,2}
Make/O Mn_LMM2 = {944,0.000,0.232,960,920,959,921,0,1,2}
Make/O Hg4d3 = {379,1.801,0.840,389,369,388,370,2,1,2}
Make/O Hg4d5 = {359,2.607,1.152,377,352,376,353,2,1,2}
Make/O Hg4f = {100,7.512,4.302,115,95,114,96,3,1,2}
Make/O Hg4f5 = {104,3.204,3.204,114,94,113,95,2,1,2}
Make/O Hg4f7 = {100,4.302,4.302,115,90,114,91,2,1,2}
Make/O Hg4p3 = {577,1.406,0.476,587,567,586,568,2,1,2}
Make/O Mo3d = {230,3.544,2.243,242,222,241,223,3,1,2}
Make/O Mo3d3 = {233,1.403,1.467,243,223,242,224,2,1,2}
Make/O Mo3d5 = {230,2.139,2.243,243,223,242,224,2,1,2}
Make/O Mo3p1 = {413,0.708,0.413,423,403,422,404,2,1,2}
Make/O Mo3p3 = {396,1.393,0.869,410,390,409,391,2,1,2}
Make/O Mo3s = {508,0.620,0.135,518,498,517,499,0,1,2}
Make/O Mo4p = {36,0.517,0.135,46,26,45,27,0,1,2}
Make/O Nd3d = {980,12.000,4.373,1025,970,1024,971,2,1,2}
Make/O Nd3d5 = {981,7.343,4.373,1000,970,999,971,3,1,2}
Make/O Nd4d = {120,2.708,2.708,133,113,132,114,2,1,2}
Make/O Ne1s = {863,1.215,1.215,876,856,875,857,3,1,2}
Make/O Ne2s = {41,0.098,0.098,55,35,54,36,0,1,2}
Make/O Ne_KLL = {669,0.000,0.461,675,655,674,656,0,1,2}
Make/O Ni2p = {855,3.702,2.468,894,844,893,845,0,1,2}
Make/O Ni2p1 = {873,1.154,1.191,883,863,882,864,0,1,2}
Make/O Ni2p3 = {855,2.309,2.468,873,848,872,849,3,1,2}
Make/O Ni3p = {67,0.693,0.416,80,60,79,61,2,1,2}
Make/O Ni3s = {112,0.200,0.136,122,102,121,103,0,1,2}
Make/O Ni_LMM = {641,0.000,1.047,660,610,659,611,2,1,2}
Make/O Ni_LMM1 = {712,0.000,0.586,725,690,724,691,0,1,2}
Make/O Ni_LMM2 = {778,0.000,0.305,785,750,784,751,0,1,2}
Make/O Nb3d = {206,3.127,2.046,217,197,216,198,3,1,2}
Make/O Nb3d3 = {209,1.249,1.314,219,199,218,200,2,1,2}
Make/O Nb3d5 = {206,1.876,2.046,217,197,216,198,2,1,2}
Make/O Nb3p1 = {379,0.660,0.432,389,369,388,370,2,1,2}
Make/O Nb3p3 = {364,1.330,0.871,377,357,376,356,2,1,2}
Make/O Nb3s = {470,0.590,0.149,480,460,479,461,0,1,2}
Make/O Nb4p = {31,0.505,0.149,41,21,40,22,0,1,2}
Make/O N1s = {402,0.499,0.499,411,391,410,392,3,1,2}
Make/O N_KLL = {1107,0.000,0.036,1090,1120,1089,1121,0,1,2}
Make/O Os4d = {279,3.653,1.134,312,272,311,273,0,1,2}
Make/O Os4d3 = {294,1.439,0.801,304,284,303,285,2,1,2}
Make/O Os4d5 = {279,2.205,1.134,312,272,311,273,2,1,2}
Make/O Os4f = {52,4.880,3.067,66,46,65,47,3,1,2}
Make/O Os4p1 = {549,0.374,0.110,559,539,558,540,2,1,2}
Make/O Os4p3 = {475,0.995,0.412,485,465,484,466,2,1,2}
Make/O O1s = {531,0.733,0.733,543,523,542,524,3,1,2}
Make/O O2s = {23,0.035,0.035,35,15,34,16,0,1,2}
Make/O O_KLL = {978,0.000,0.171,990,960,989,961,2,1,2}
Make/O Pd3d = {337,5.637,3.309,350,330,349,331,3,1,2}
Make/O Pd3d3 = {342,2.323,2.323,352,332,351,333,2,1,2}
Make/O Pd3d5 = {337,3.309,3.309,352,332,351,333,2,1,2}
Make/O Pd3p1 = {561,0.787,0.404,571,551,570,552,2,1,2}
Make/O Pd3p3 = {534,1.586,0.821,548,528,547,529,2,1,2}
Make/O Pd3s = {671,0.740,0.821,681,661,680,662,0,1,2}
Make/O Pd4p = {52,0.570,0.821,62,42,61,43,0,1,2}
Make/O Pd_MNN = {1159,0.000,0.942,1170,1140,1169,1141,0,1,2}
Make/O P2p = {133,0.525,0.485,143,123,142,124,3,1,2}
Make/O P2s = {191,0.380,0.276,205,185,204,186,2,1,2}
Make/O Pt4d3 = {333,1.596,0.614,343,323,342,324,2,1,2}
Make/O Pt4d5 = {316,2.377,1.114,333,308,332,309,2,1,2}
Make/O Pt4f = {73,6.080,3.524,86,66,85,67,3,1,2}
Make/O Pt4f5 = {76,2.553,2.553,86,66,85,67,2,1,2}
Make/O Pt4f7 = {73,3.524,3.524,88,68,87,69,2,1,2}
Make/O Pt4p3 = {521,1.083,0.413,531,511,530,512,2,1,2}
Make/O Pt_NNN = {1317,0.000,0.000,1350,1300,1349,1301,0,1,2}
Make/O K2p = {293,1.552,1.039,307,287,306,288,3,1,2}
Make/O K2p1 = {296,0.512,0.512,306,286,305,287,2,1,2}
Make/O K2p3 = {293,1.039,1.039,307,287,306,288,2,1,2}
Make/O K2s = {378,0.517,0.312,392,372,391,373,0,1,2}
Make/O K_LMM = {1239,0.000,0.304,1250,1230,1249,1231,0,1,2}
Make/O Pr3d = {930,12.000,6.335,970,920,969,921,2,1,2}
Make/O Pr3d5 = {932,7.343,6.335,957,927,956,928,3,1,2}
Make/O Pr4d = {114,2.715,2.715,127,107,126,108,2,1,2}
Make/O Re4d = {263,3.491,1.142,295,255,294,256,0,1,2}
Make/O Re4d3 = {277,1.374,0.794,287,267,286,268,2,1,2}
Make/O Re4d5 = {263,2.107,1.142,295,255,294,256,2,1,2}
Make/O Re4f = {43,4.339,2.939,55,35,54,36,3,1,2}
Make/O Re4p1 = {521,0.368,0.123,531,511,530,512,2,1,2}
Make/O Re4p3 = {449,0.962,0.429,459,439,458,440,2,1,2}
Make/O Re4s = {628,0.470,0.094,638,618,637,619,0,1,2}
Make/O Rh3d = {309,5.092,2.981,322,302,321,303,3,1,2}
Make/O Rh3d3 = {314,2.106,2.106,324,304,323,305,2,1,2}
Make/O Rh3d5 = {309,2.816,2.816,324,304,323,305,2,1,2}
Make/O Rh3p1 = {522,0.770,0.413,532,512,531,513,2,1,2}
Make/O Rh3p3 = {498,1.560,0.836,513,493,512,494,2,1,2}
Make/O Rh3s = {629,0.710,0.836,639,619,638,620,0,1,2}
Make/O Rh4p = {48,0.570,0.836,58,38,57,39,0,1,2}
Make/O Rh_MNN = {1185,0.000,0.724,1205,1170,1204,1171,0,1,2}
Make/O Rb3d = {110,1.673,1.455,125,105,124,106,3,1,2}
Make/O Rb3p = {238,1.669,0.860,260,235,259,236,0,1,2}
Make/O Rb3p1 = {247,0.550,0.422,257,237,256,238,2,1,2}
Make/O Rb3p3 = {238,1.117,0.860,260,235,259,236,2,1,2}
Make/O Rb3s = {322,0.470,0.222,332,312,331,313,0,1,2}
Make/O Ru3d = {282,4.529,2.705,295,275,294,276,3,1,2}
Make/O Ru3d3 = {286,1.821,1.821,296,276,295,277,2,1,2}
Make/O Ru3d5 = {282,2.705,2.705,297,277,296,278,2,1,2}
Make/O Ru3p1 = {485,0.740,0.410,495,475,494,476,2,1,2}
Make/O Ru3p3 = {463,1.496,0.840,478,458,477,459,2,1,2}
Make/O Ru3s = {586,0.680,0.840,596,576,595,577,0,1,2}
Make/O Ru4p = {43,0.550,0.840,53,33,52,34,0,1,2}
Make/O Ru_MNN = {1212,0.000,0.510,1230,1200,1229,1201,0,1,2}
Make/O Sm3d3 = {1108,5.001,2.694,1130,1091,1129,1092,2,1,2}
Make/O Sm3d5 = {1083,7.343,2.309,1100,1065,1099,1066,3,1,2}
Make/O Sm4d = {132,2.694,2.694,147,127,146,128,2,1,2}
Make/O Sc2p = {402,1.959,1.306,414,394,413,395,3,1,2}
Make/O Sc2p1 = {407,0.651,0.497,417,397,416,398,2,1,2}
Make/O Sc2p3 = {402,1.306,1.306,417,397,416,398,2,1,2}
Make/O Sc2s = {501,0.564,0.293,516,496,515,497,0,1,2}
Make/O Sc_LMM = {1149,0.000,0.175,1160,1135,1159,1136,0,1,2}
Make/O Sc_LMM1 = {1118,0.000,0.375,1128,1100,1127,1101,0,1,2}
Make/O Se3d = {56,0.821,0.640,69,49,68,50,3,1,2}
Make/O Se3p1 = {169,0.498,0.064,179,159,178,160,0,1,2}
Make/O Se3p3 = {163,0.996,0.296,173,153,172,154,0,1,2}
Make/O Se3s = {232,0.380,0.000,242,222,241,223,0,1,2}
Make/O Se_LMM = {300,0.000,1.128,182,157,181,158,2,1,2}
Make/O Se_LMM1 = {184,0.000,0.500,202,177,201,178,0,1,2}
Make/O Se_LMM2 = {143,0.000,2.355,163,133,162,134,0,1,2}
Make/O Se_LMM3 = {259,0.000,0.322,279,249,278,250,0,1,2}
Make/O Se_LMM4 = {288,0.000,0.397,308,278,307,279,0,1,2}
Make/O Si2p = {102,0.368,0.341,114,94,113,95,3,1,2}
Make/O Si2s = {153,0.347,0.254,165,145,164,146,2,1,2}
Make/O Ag3d = {368,6.277,3.742,382,362,381,363,3,1,2}
Make/O Ag3d3 = {374,2.527,2.527,374,354,373,355,2,1,2}
Make/O Ag3d5 = {368,3.742,3.742,384,364,383,365,2,1,2}
Make/O Ag3p1 = {604,0.796,0.435,614,594,613,595,2,1,2}
Make/O Ag3p3 = {573,1.649,0.814,585,565,584,566,2,1,2}
Make/O Ag3s = {719,0.770,0.814,729,709,728,710,2,1,2}
Make/O Ag4d = {5,0.664,2.422,15,-5,14,-4,2,1,2}
Make/O Ag4p = {60,0.580,0.814,70,50,69,51,2,1,2}
Make/O Ag4s = {98,0.178,0.814,108,88,107,89,2,1,2}
Make/O Ag_MNN = {1129,0.000,0.750,1145,1120,1144,1121,2,1,2}
Make/O Ag_MNN1 = {1135,0.000,0.946,1145,1120,1144,1121,0,1,2}
Make/O Na1s = {1072,1.102,1.102,1086,1066,1085,1067,3,1,2}
Make/O Na2p = {31,0.084,0.070,41,21,40,22,0,1,2}
Make/O Na2s = {64,0.180,0.166,78,58,77,59,0,1,2}
Make/O Na_KLL = {493,0.000,1.440,520,488,519,489,2,1,2}
Make/O Na_KLL1 = {532,0.000,0.330,545,528,544,529,0,1,2}
Make/O Sr3d = {133,1.992,1.669,153,128,152,129,3,1,2}
Make/O Sr3p = {269,1.747,0.873,291,266,290,267,0,1,2}
Make/O Sr3p1 = {280,0.579,0.416,290,270,289,271,2,1,2}
Make/O Sr3p3 = {269,1.164,0.873,291,266,290,267,2,1,2}
Make/O Sr3s = {358,0.500,0.146,368,348,367,349,0,1,2}
Make/O Sr4p = {21,0.461,0.146,31,11,30,12,0,1,2}
Make/O S2p = {165,0.717,0.650,175,155,174,156,3,1,2}
Make/O S2s = {229,0.426,0.309,242,222,241,223,2,1,2}
Make/O S_LMM = {1336,0.000,0.000,1350,1320,1349,1321,0,1,2}
Make/O Ta4d = {229,3.224,1.160,260,220,259,221,0,1,2}
Make/O Ta4d3 = {241,1.282,0.769,251,231,250,232,2,1,2}
Make/O Ta4d5 = {229,1.934,1.160,260,220,259,221,2,1,2}
Make/O Ta4f = {25,3.384,2.665,37,17,36,16,3,1,2}
Make/O Ta4p1 = {464,0.391,0.161,474,454,473,455,2,1,2}
Make/O Ta4p3 = {403,0.820,0.475,413,393,412,394,2,1,2}
Make/O Ta4s = {566,0.440,0.131,576,556,575,557,0,1,2}
Make/O Te3d3 = {585,4.088,4.088,595,575,594,576,2,1,2}
Make/O Te3d5 = {575,5.850,5.850,587,567,586,568,3,1,2}
Make/O Te3p1 = {873,0.794,0.397,883,863,882,864,2,1,2}
Make/O Te3p3 = {822,1.750,0.875,832,812,831,813,2,1,2}
Make/O Te4d = {43,1.721,1.358,57,37,56,38,2,1,2}
Make/O Te_MNN = {995,0.000,0.883,1015,985,1014,986,2,1,2}
Make/O Te_MNN1 = {1005,0.000,1.011,1015,985,1014,986,0,1,2}
Make/O Tb3d = {1244,12.000,0.621,1275,1235,1274,1236,2,1,2}
Make/O Tb4d = {150,2.672,2.672,170,140,169,141,3,1,2}
Make/O Tl4d3 = {406,1.765,0.782,416,396,415,395,2,1,2}
Make/O Tl4d5 = {385,2.693,1.137,403,378,402,379,2,1,2}
Make/O Tl4f = {118,8.334,4.743,132,112,131,113,3,1,2}
Make/O Tl4f5 = {122,3.585,3.585,132,112,131,113,2,1,2}
Make/O Tl4f7 = {118,4.743,4.743,133,113,132,114,2,1,2}
Make/O Tl4p3 = {608,1.164,0.434,618,598,617,599,2,1,2}
Make/O Tl5d = {13,1.276,0.780,28,8,27,9,0,1,2}
Make/O Th4d3 = {714,2.252,0.783,724,704,723,705,2,1,2}
Make/O Th4d5 = {677,3.524,1.208,695,670,694,671,2,1,2}
Make/O Th4f5 = {344,7.450,7.450,354,334,353,335,2,1,2}
Make/O Th4f7 = {335,9.568,9.568,347,327,346,328,3,1,2}
Make/O Th5d = {87,2.060,1.236,102,82,101,83,0,1,2}
Make/O Th5d3 = {94,0.822,0.822,104,84,103,85,0,1,2}
Make/O Th5d5 = {87,1.236,1.236,102,82,101,83,0,1,2}
Make/O Tm4d = {180,2.637,2.637,185,165,184,166,3,1,2}
Make/O Sn3d = {488.5,8.288,4.890,500,480,499,481,2,1,2}
Make/O Sn3d3 = {494,3.398,3.398,504,484,503,485,2,1,2}
Make/O Sn3d5 = {486,4.890,4.890,500,480,499,481,3,1,2}
Make/O Sn3p1 = {757,0.824,0.350,767,747,766,748,2,1,2}
Make/O Sn3p3 = {715,1.732,0.734,730,710,729,711,2,1,2}
Make/O Sn3s = {885,0.870,0.804,895,875,894,876,0,1,2}
Make/O Sn4d = {25,1.100,0.804,35,15,34,16,2,1,2}
Make/O Sn_MNN = {1049,0.000,0.814,1070,1040,1069,1041,2,1,2}
Make/O Sn_MNN1 = {1058,0.000,0.902,1070,1040,1069,1041,0,1,2}
Make/O Ti2p = {458,2.077,1.385,469,449,468,450,3,1,2}
Make/O Ti2p1 = {464,0.690,0.690,474,454,473,455,2,1,2}
Make/O Ti2p3 = {458,1.385,1.385,476,451,475,452,2,1,2}
Make/O Ti2s = {565,0.589,0.284,578,558,577,559,2,1,2}
Make/O Ti3p = {37,0.090,0.210,52,32,51,33,2,1,2}
Make/O Ti3s = {59,0.030,0.171,70,50,69,51,0,1,2}
Make/O Ti_LMM = {1068,0.000,0.332,1080,1055,1079,1056,0,1,2}
Make/O Ti_LMM1 = {1098,0.000,0.182,1110,1080,1109,1081,0,1,2}
Make/O W4d = {245,3.328,1.152,276,236,275,237,0,1,2}
Make/O W4d3 = {257,1.387,0.827,267,247,266,248,2,1,2}
Make/O W4d5 = {245,2.010,1.152,276,236,275,237,2,1,2}
Make/O W4f = {34,3.863,2.810,46,26,45,27,3,1,2}
Make/O W4p1 = {491,0.386,0.125,501,481,500,482,2,1,2}
Make/O W4p3 = {425,0.830,0.469,435,415,434,416,2,1,2}
Make/O W4s = {594,0.460,0.118,604,584,603,585,0,1,2}
Make/O U4d3 = {781,2.359,0.836,791,771,790,772,2,1,2}
Make/O U4d5 = {739,3.697,1.248,757,732,756,733,2,1,2}
Make/O U4f5 = {391,8.474,8.474,401,381,400,382,2,1,2}
Make/O U4f7 = {380,10.800,10.800,387,372,386,373,3,1,2}
Make/O U5d = {96,2.189,1.368,110,90,109,91,0,1,2}
Make/O U5d3 = {104,0.818,0.818,114,94,113,95,0,1,2}
Make/O U5d5 = {96,1.368,1.368,110,90,109,91,0,1,2}
Make/O V2p = {515,2.184,1.456,536,511,535,512,0,1,2}
Make/O V2p1 = {523,0.725,0.725,533,513,532,514,0,1,2}
Make/O V2p3 = {515,1.456,1.456,527,507,526,508,3,1,2}
Make/O V2s = {630,0.688,0.134,640,620,639,621,0,1,2}
Make/O V3p = {40,0.140,0.224,53,33,52,34,0,1,2}
Make/O V3s = {66,0.050,0.183,75,55,74,56,0,1,2}
Make/O V_LMM = {1014,0.000,0.114,1030,1000,1029,1001,0,1,2}
Make/O V_LMM1 = {1048,0.000,0.413,1060,1040,1059,1041,0,1,2}
Make/O V_LMM2 = {977,0.000,0.307,1000,955,999,956,0,1,2}
Make/O Xe3d3 = {685,4.404,4.404,695,675,694,676,2,1,2}
Make/O Xe3d5 = {672,6.670,6.670,684,664,683,665,3,1,2}
Make/O Xe4d = {63,2.387,1.610,77,57,76,58,2,1,2}
Make/O Xe_MNN = {955,0.000,1.127,970,920,969,921,2,1,2}
Make/O Yb4d = {185,2.631,2.631,215,165,214,166,3,1,2}
Make/O Yb4p1 = {389,0.385,0.734,405,379,404,380,2,1,2}
Make/O Yb4p3 = {333,0.760,0.734,348,328,347,329,2,1,2}
Make/O Yb4s = {482,0.410,0.734,492,472,491,473,0,1,2}
Make/O Y3d = {158,2.343,1.824,171,151,170,152,3,1,2}
Make/O Y3p = {301,1.833,0.885,323,298,322,299,0,1,2}
Make/O Y3p1 = {313,0.610,0.434,323,303,322,304,2,1,2}
Make/O Y3p3 = {301,1.222,0.885,323,298,322,299,2,1,2}
Make/O Y3s = {395,0.530,0.143,405,385,404,386,0,1,2}
Make/O Y4p = {24,0.460,0.143,34,14,33,15,0,1,2}
Make/O Zn2p1 = {1045,1.384,1.253,1055,1035,1054,1036,2,1,2}
Make/O Zn2p3 = {1022,2.768,2.680,1036,1016,1035,1017,3,1,2}
Make/O Zn3d = {10,0.371,0.549,20,0,19,1,0,1,2}
Make/O Zn3p = {89,1.029,0.549,104,84,103,85,2,1,2}
Make/O Zn3s = {140,0.260,0.549,150,130,149,131,0,1,2}
Make/O Zn_LMM = {495,2.911,0.546,515,480,514,481,2,1,2}
Make/O Zn_LMM1 = {472,0.000,1.919,492,480,491,481,0,1,2}
Make/O Zn_LMM2 = {582,0.000,0.242,600,570,599,571,0,1,2}
Make/O Zn_LMM3 = {573,0.000,0.379,600,570,599,571,0,1,2}
SetDataFolder savDF
end

Function VersaProbeConcentration(peakAreas, elementList)
Wave peakAreas
String elementList

Make/O/N=(numpnts(peakAreas)) corrAreas, percentComp
Variable i

For(i=0;i<itemsInList(elementList);i+=1)
Wave w1 = $("root:SensitivityFactors:" + stringFromList(i, elementList))

if(!WaveExists($("root:crossSections:" + stringFromList(i, elementList))))
	getCrossSections(w1)
endif

corrAreas[i] = peakAreas[i]/w1[1]
endfor

For(i=0;i<itemsInList(elementList);i+=1)
	Wave w2 = $("root:crossSections:" + stringFromList(i, elementList))
	//corrAreas[i] = (corrAreas[i]*(4*pi*(1-1/4*w2[4]*(3*cos(pi/180*45)^2-1))))
endfor
	percentComp = corrAreas/sum(corrAreas)*100
end

Function batchVersaProbeConcentration(peakAreasWaveList,elementList,ignoreCarbon)
String peakAreasWaveList,elementList
Variable ignoreCarbon
Variable p1
NewDataFolder/O root:Batch_VersaProbe_Comp

For(p1 = 0; p1 < itemsInList(peakAreasWaveList); p1 += 1)
	Wave tempWave = $stringFromList(p1,peakAreasWaveList)
	NewDataFolder/O root:Batch_VersaProbe_Comp:$NameOfWave(tempWave)
	String newAreaWave = "root:Batch_VersaProbe_Comp:" + NameOfWave(tempWave) + ":areas"
	Make/O/N=(numpnts(tempWave)) $newAreaWave
	Wave AreasWave = $newAreaWave
	AreasWave = tempWave
	if(ignoreCarbon)
		AreasWave[0] = 0
	endif
endFor
end

Function angCorrVersaProbe()
String savDF = getDataFolder(1)
NewDataFolder/O/S root:VersaProbe_Angular_Correction
Make/O/N=5000 angCorr
SetScale/I x,0,90,angCorr
angCorr = -57718 + 7760.4*x + -249.96*x^2 + 5.1771*x^3 + -0.053105*x^4 + 0.00020873*x^5
Variable maxVal = WaveMax(angCorr)
angCorr = angCorr/maxVal
SetDataFolder savDF
end

Function xpsCorrAngVP(inputWave)
Wave inputWave
Wave angCorr = root:VersaProbe_Angular_Correction:angCorr
String waveNote = note(inputWave)
Variable corrValue = str2num(stringFromList(4,StringByKey("StagePosition", waveNote, ":", "\r")," "))
Duplicate/O inputWave, $(GetWavesDataFolder(inputWave,2) + "_angCorr")
Wave w1 = $(GetWavesDataFolder(inputWave,2) + "_angCorr")
w1 = w1/angCorr(corrValue)
end

Function xpsVPalign(inputWaveList, peakLoc, startPnt, endPnt)
String inputWaveList
Variable peakLoc, startPnt, endPnt

Variable numWavesInList = ItemsInList(inputWaveList)
Variable lowX = inf, highX = inf, p1

Do
	Wave tempWave = $(StringFromList(p1, inputWaveList))
	WaveStats/Q/R=(pnt2x(tempWave,endPnt), pnt2x(tempWave,startPnt)) tempWave
	Variable shiftValue = V_maxLoc - peakLoc
	Duplicate/O tempWave, $(GetWavesDataFolder(tempWave,2)  + "_faCorr")
	Wave w1 = $(GetWavesDataFolder(tempWave,2)  + "_faCorr")
	Setscale/I x,(pnt2x(tempWave,0) - shiftValue),(pnt2x(tempWave,numpnts(tempWave)-1) - shiftValue), w1
	p1+=1
while(p1 < numWavesInList)

p1 = 0
Do
	Wave tempWave = $(StringFromList(p1, inputWaveList) + "_faCorr")
	if(lowX == inf)
		lowX = pnt2x(tempWave, numpnts(tempWave)-1)
	endif
	if(lowX < pnt2x(tempWave, numpnts(tempWave)-1))
		lowX = pnt2x(tempWave, numpnts(tempWave)-1)
	endif
	p1+=1
while(p1 < numWavesInList)

print lowX
p1=0
Do
	Wave tempWave = $(StringFromList(p1, inputWaveList) + "_faCorr")
	Print GetWavesDataFolder(tempWave,2), x2pnt(tempWave,lowX)
	if((numpnts(tempWave) -1) - x2pnt(tempWave,lowX) > 0)
		DeletePoints x2pnt(tempWave,lowX),(numpnts(tempWave) - 1) - x2pnt(tempWave,lowX) , tempWave
	endif

	p1+=1
while(p1 < numWavesInList)

p1=0
Do
	Wave tempWave = $(StringFromList(p1, inputWaveList) + "_faCorr")
	if(highX == inf)
		highX = pnt2x(tempWave, 0)
	endif
	if(highX > pnt2x(tempWave, 0))
		highX = pnt2x(tempWave, 0)
	endif
	p1+=1
while(p1 < numWavesInList)
Print highX

p1=0
Do
	Wave tempWave = $(StringFromList(p1, inputWaveList) + "_faCorr")
	Print GetWavesDataFolder(tempWave,2),x2pnt(tempWave,highX)
	if(x2pnt(tempWave,highX) > 0)
		DeletePoints 0,x2pnt(tempWave,highX), tempWave
	endif
	Setscale/I x,highX,lowX, tempWave
	p1+=1
while(p1 < numWavesInList)

end

//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//
//																	 //
//																	 //
//  				        	Normalize Wave List							 //
//																	 //
//																	 //
//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//

Function normWave(inputWave)
Wave inputWave
String wName = NameOfWave(inputWave)

Make/O/N=(numpnts(inputWave)) $(wName + "_N")
Wave w1 = $(wName + "_N")

w1 = (inputWave - WaveMin(inputWave))/(WaveMax(inputWave) - Wavemin(inputWave))

end

Function interWaves(inputWave)
Wave inputWave
String wName = NameOfWave(inputWave)
Variable/G $(wName + "_Sum")
NVAR intRes = $(wName + "_Sum")

Make/O/N=(numpnts(inputWave)) $(wName +"_Int")
 Wave w1 = $(wName +"_Int")
 w1 = inputWave
 Variable startPnt = pnt2x(w1,0), endPnt = pnt2x(w1,30)
WaveStats/Q/R=(startPnt,endPnt) w1

w1 = w1- V_avg

intRes = sum(w1)
end

Function roughComp(elementList)
String elementList
Variable/G total = 0
Variable p1

Do
NVAR temp = $(stringfromlist(p1, elementList))
total += temp
p1 += 1
while(p1 < itemsinlist(elementList))

p1=0

Do
Variable/G $(stringfromlist(p1, elementList) + "_percent")
NVAR temp = $(stringfromlist(p1, elementList) + "_percent")
NVAR temp2 = $(stringfromlist(p1, elementList))
temp = (temp2/total) * 100
p1 += 1
while(p1 < itemsinlist(elementList))

end

Function reportRoughComp(elementList)
String elementList
String TopWin = UniqueName("Composition", 10, 1)
NewNoteBook/F=0/K=1/N=$TopWin

NoteBook kwTopWin, text = TopWin + "\r"
Variable p1 = 0
Do
NVAR temp =  $(stringfromlist(p1, elementList) + "_percent")
NoteBook kwTopWin, text= stringfromlist(p1, elementList)[0,strsearch(stringfromlist(p1, elementList), "_", 0) -1] + ": " + num2str(temp) + "%\r"

Make/O/N=18 $("root:" + stringfromlist(p1, elementList)[0,strsearch(stringfromlist(p1, elementList), "_", 0) -1])
Wave w1 = $("root:" + stringfromlist(p1, elementList)[0,strsearch(stringfromlist(p1, elementList), "_", 0) -1])

w1[str2num(TopWin[strsearch(TopWin, "n", 0)+1,inf])-1] = temp

p1+=1
while(p1 < itemsinlist(elementList))

end

Function WalkDataFoldersDisp(path)
	String path				// to walk all data folders, pass "root:" WITH the : char. To use the current data folder, use ":" or GetDataFolder(1)

	String savDF= GetDataFolder(1)
	SetDataFolder path
	
	Wave w1 = $("Ga2p_subtracted")
	Variable test = waveExists(w1)
	If(test)
		Wave w1 = $"C1s_subtracted"
		Wave w2 = $"Ga2p_subtracted"
		Wave w3 = $"In3d5_subtracted"
		Wave w4 = $"O1s_subtracted"
		Wave w5 = $"Zn2p3_subtracted"
		Display w1,w2,w3,w4,w5
	endif

	Variable i, numDataFolders = CountObjects(":",4)
	for(i=0; i<numDataFolders; i+=1)
		String nextPath = GetIndexedObjName(":",4,i)
		WalkDataFoldersDisp(nextPath)
	endfor

	SetDataFolder savDF
End
Function linExp(w,x) : FitFunc
	Wave w
	Variable x

	//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
	//CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
	//CurveFitDialog/ Equation:
	//CurveFitDialog/ f(x) = C1*(C2*x+C3)^C4*exp(-(C2*x-C3)/C5)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 5
	//CurveFitDialog/ w[0] = C1
	//CurveFitDialog/ w[1] = C2
	//CurveFitDialog/ w[2] = C3
	//CurveFitDialog/ w[3] = C4
	//CurveFitDialog/ w[4] = C5

	return w[0]*(w[1]*x+w[2])^w[3]*exp(-(w[1]*x-w[2])/w[4])
End

Function WalkDataFoldersDup(path, folderPath, waveListMatchStr, waveListOptionsStr)
	String path				// to walk all data folders, pass "root:" WITH the : char. To use the current data folder, use ":" or GetDataFolder(1)
	String folderPath
	String waveListMatchStr		// "*" or "wave*" or "!*x", etc. See WaveList for details.
	String waveListOptionsStr	// usually "". see WaveList's optionsStr for details.
	
	String savDF= GetDataFolder(1)
	SetDataFolder path
	
	NewDataFolder/O $("root:" + folderPath)
	
	String listOfMatchingWaves= WaveList(waveListMatchStr,";",waveListOptionsStr)
	Variable j = 0
	for(j = 0; j < itemsInList(listOfMatchingWaves, ";"); j+=1)
		String tempDF = GetDataFolder(1)
		String tempName
		String tempFolderName = GetDataFolder(0)
		SetDataFolder $("root:" + folderPath)
		//tempFolderName = tempFolderName[0,strsearch(tempFolderName,"_",inf,1)-1]
		//tempFolderName = tempFolderName[0,strsearch(tempFolderName,"_",inf,1)-1]
		
		//tempName = tempFolderName + "_" + stringFromList(j, listOfMatchingWaves, ";")
		//tempName = stringFromList(j, listOfMatchingWaves, ";")
		//tempName = replaceString("_sub_smth",tempName,"_ss")
		//tempName = UniqueName(replaceString("GIZO_", tempName, "") + "_",1,1)
		tempName = UniqueName(stringFromList(j, listOfMatchingWaves, ";") + "_",1,1)
		SetDataFolder tempDF
		Duplicate $(stringFromList(j, listOfMatchingWaves, ";")) $("root:" + folderPath +":"+ tempName)
	endfor

	Variable i, numDataFolders = CountObjects(":",4)
	for(i=0; i<numDataFolders; i+=1)
		String nextPath = GetIndexedObjName(":",4,i)
		WalkDataFoldersDup(nextPath, folderPath,waveListMatchStr, waveListOptionsStr)
	endfor

	SetDataFolder savDF
End

function replaceStrWaves(inputWave, stringRemove, stringReplace)
Wave inputWave
String stringRemove,stringReplace

String oldName = NameOfWave(inputWave)
String newName = ReplaceString(stringRemove, oldName, stringReplace)
Rename $oldName, $newName
	
end

function changeName(waveListMatchStr,waveListOptionsStr)
	String waveListMatchStr,waveListOptionsStr
	String tempName
	String listOfMatchingWaves= WaveList(waveListMatchStr,";",waveListOptionsStr) 
	
	Variable j = 0
	for(j = 0; j < itemsInList(listOfMatchingWaves, ";"); j+=1)
		tempName = stringFromList(j, listOfMatchingWaves, ";")
		Variable num = str2num(tempName[strsearch(tempName,"_",inf,1)+1,inf])

		if(mod(num,3) == 1)
			tempName = stringFromList(j, listOfMatchingWaves, ";") + "_45deg"
			Duplicate/O $(stringFromList(j, listOfMatchingWaves, ";")) $(tempName)
			Wave w1 = $tempName
			Wave w2 = root:'20130117_Au':Au_ang:angScaleVB
			w1 = w1/w2(45)
		elseif(mod(num,3) == 2)
			tempName = stringFromList(j, listOfMatchingWaves, ";") + "_30deg"
			Duplicate/O $(stringFromList(j, listOfMatchingWaves, ";")) $(tempName)
			Wave w1 = $tempName
			Wave w2 = root:'20130117_Au':Au_ang:angScaleVB
			w1 = w1/w2(30)
		elseif(mod(num,3) == 0)
			tempName = stringFromList(j, listOfMatchingWaves, ";") + "_80deg"
			Duplicate/O $(stringFromList(j, listOfMatchingWaves, ";")) $(tempName)
			Wave w1 = $tempName
			Wave w2 = root:'20130117_Au':Au_ang:angScaleVB
			w1 = w1/w2(80)
		endif
		
	endfor
end

Function angVar(w,theta) : FitFunc
	Wave w
	Variable theta

	//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
	//CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
	//CurveFitDialog/ Equation:
	//CurveFitDialog/ f(theta) = I0*(exp(-(d/(lamda*sin(theta)))))
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ theta
	//CurveFitDialog/ Coefficients 3
	//CurveFitDialog/ w[0] = d
	//CurveFitDialog/ w[1] = lamda
	//CurveFitDialog/ w[2] = I0

	return w[2]*(exp(-(w[0]/(w[1]*sin(theta)))))
End

function depthIMFP(lowBindWaveClean, highBindWaveClean, lowBindWaveDirty, highBindWaveDirty, angleDeg)
Wave lowBindWaveClean, highBindWaveClean, lowBindWaveDirty, highBindWaveDirty
Variable angleDeg

Variable angleRad = cos(angleDeg*3.14159265358979323846264338327950288419716939937510/180)
Variable highLamda, lowLamda

Print  recursiveStrSearch(NameOfWave(lowBindWaveClean), "_", 4, 0, inf, -1)
Print "Angle: ", angleDeg, "degrees"
Wavestats/Q lowBindWaveClean
lowLamda = (49/(1486.7 - V_maxloc)^2+0.11*(1486.7 - V_maxloc)^.5)*13.5
Print "Low BE Peak Location: ", V_maxloc, "eV"
Wavestats/Q highBindWaveClean
highLamda = (49/(1486.7 - V_maxloc)^2+0.11*(1486.7 - V_maxloc)^.5)*13.5
Print "High BE Peak Location: ", V_maxloc, "eV"

Print "Low Lamda: ", lowLamda, "Angstroms"
Print "High Lamda: ", highLamda, "Angstroms"

Variable thickness = (lowLamda*highLamda)/(highLamda-lowLamda)*sin(angleRad)*ln((Wavemax(lowBindWaveClean)/Wavemax(highBindWaveClean))/((Wavemax(lowBindWaveDirty)/Wavemax(highBindWaveDirty))))

Print "Overlayer Thickness: ", thickness, "Angstroms"

end

Function/S recursiveStrSearch(inputStr, searchStr, numOfTimes, stringStart, searchStart, direction)
String inputStr,searchStr
Variable stringStart, numOfTimes,  searchStart,  direction
Variable i
for(i=0;i<numOfTimes;i+=1)
	 inputStr = inputStr[stringStart,strSearch(inputStr,searchStr,searchStart,direction)-1]
endfor

Return inputStr
end

Function WalkDataFoldersDup2(path, folderPath, waveListMatchStr, waveListOptionsStr)
	String path				// to walk all data folders, pass "root:" WITH the : char. To use the current data folder, use ":" or GetDataFolder(1)
	String folderPath
	String waveListMatchStr		// "*" or "wave*" or "!*x", etc. See WaveList for details.
	String waveListOptionsStr	// usually "". see WaveList's optionsStr for details.
	
	String savDF= GetDataFolder(1)
	NewDataFolder/O $("root:" + folderPath)
	SetDataFolder path
	
	String listOfMatchingWaves= WaveList(waveListMatchStr,";",waveListOptionsStr)
	Variable j = 0
	if(StringMatch(path,"factors*"))
		for(j = 0; j < itemsInList(listOfMatchingWaves, ";"); j+=1)
			Duplicate $(stringFromList(j, listOfMatchingWaves, ";")) $("root:" + folderPath +":"+ stringFromList(j, listOfMatchingWaves, ";"))
		endfor
	endif

	Variable i, numDataFolders = CountObjects(":",4)
	for(i=0; i<numDataFolders; i+=1)
		String nextPath = GetIndexedObjName(":",4,i)
		WalkDataFoldersDup2(nextPath, folderPath,waveListMatchStr, waveListOptionsStr)
	endfor

	SetDataFolder savDF
End

function changeName2(waveListMatchStr,waveListOptionsStr)
	String waveListMatchStr,waveListOptionsStr
	String tempName
	String listOfMatchingWaves= WaveList(waveListMatchStr,";",waveListOptionsStr) 
	
	Variable j = 0
	for(j = 0; j < itemsInList(listOfMatchingWaves, ";"); j+=1)
		tempName = ReplaceString("VBXPSangCorr_",stringfromlist(j,listofMatchingWaves,";"), "")
		//tempName = ReplaceString("Waves",tempName, "")
		Rename $stringfromlist(j,listofMatchingWaves,";"), $tempName
	endfor
end

Function ScaleAuFunct()
Variable numWaves = CountObjects(":",1)
Variable i = 0
Wave w2 =  root:'20130117_Au':Au_ang:angScale
for(i = 0; i < numWaves; i+=1)
	String nameWave = GetIndexedObjName(":",1,i)
	Wave w = $nameWave
	Variable angle = (i*2+30)
	w = w / w2(angle)
endfor
end

Function ManualShift(inputWaveList,inputShiftWave)
String inputWaveList
Wave inputShiftWave
Variable p1 = 0

Do
	Wave w1 = $stringFromList(p1, inputWaveList)
	Duplicate w1, $(GetWavesDataFolder(w1,2) + "_shift")
	Wave w2 = $(GetWavesDataFolder(w1,2) + "_shift")
	SetScale x,pnt2x(w2,0)+inputShiftWave[p1],pnt2x(w2, numpnts(w2))+inputShiftWave[p1],w2
	p1+=1
while(p1 < itemsInList(inputWaveList))
end

Function WaveShave(inputWaveList)
String InputWaveList 
//List of waves to work on. If you put all the waves you want to work on in one folder with no other 
//waves you can use waveShave(WaveList("*",";","")) from the command line to work on all of those waves.

Variable p1, lowX = inf, highX = inf, numWavesInList = ItemsInList(inputWaveList) 
//p1 for doing the iteration, lowX/highX to hold the x values (set to inf to make testing easier)
//numWavesInList holds the number of waves in inputWaveList

//All these Do/while loops should actually be For loops instead. This is sloppy.

Do //Makes new waves to hold the results and increases the number of points
	Wave tempWave = $(StringFromList(p1, inputWaveList)) //tempWave = an input wave from the list
	//Duplicate/O tempWave, $(GetWavesDataFolder(tempWave,2)  + "_faCorr") //Make a copy of tempWave in the same folder as tempWave with _faCorr at the end
	Duplicate/O tempWave, $(NameOfWave(tempWave)  + "_faCorr")
	Wave tempWave2 = $(StringFromList(p1, inputWaveList) + "_faCorr") //tempWave = copy of input wave from list
	Interpolate2/N=5000/T=1/Y=tempWave2/F=0 tempWave
	p1+=1
while(p1 < numWavesInList)

//Reset p1 again. Should have used For loops
p1=0
//This loop and the lowX loop could be combined into one loop. Might run slightly faster that way. This is inefficient coding.
Do  //Finds the smallest X value at the start of the waves in the list by replacing highX whenever the X value at the start of a wave is smaller 
	Wave tempWave = $(StringFromList(p1, inputWaveList) + "_faCorr")
	if(lowX == inf) //Same as before
		lowX = pnt2x(tempWave, (numpnts(tempWave)-1))
	endif
	if(lowX > pnt2x(tempWave, (numpnts(tempWave)-1))) //Same as before
		lowX = pnt2x(tempWave,(numpnts(tempWave)-1))
	endif
	p1+=1
while(p1 < numWavesInList)
//Reset p1 again. Should have used For loops
p1 = 0
//This loop could probably be combined with the loop that chops the end of the wave (as long as the end chopping happens first). Might introduce strange bugs.
Do  //Chops points off the start of the wave remembering that the second argument is the number of points and that point counting starts at zero.
     //Therefore the point of the highX is the number of points you want to delete.
	Wave tempWave = $(StringFromList(p1, inputWaveList) + "_faCorr")
	if(numpnts(tempWave) - x2pnt(tempWave,lowX) > 0)
		DeletePoints (x2pnt(tempWave,lowX)-1),(numpnts(tempWave)-x2pnt(tempWave,lowX)), tempWave
	endif
	p1+=1
while(p1 < numWavesInList)

//Reset p1 so I can use it again. Wouldn't need to do this if I used For loops. Again this is sloopy (easy to forget this and have a bug).
p1 = 0
Do //Finds the largest X value at the end of the waves in the list by replacing lowX whenever the X value at the end of a wave is bigger 
	Wave tempWave = $(StringFromList(p1, inputWaveList) + "_faCorr")
	if(highX == inf) //This is why lowX/highX were set to inf. Makes testing for the first time through the loop slightly easier
		highX = pnt2x(tempWave, 0)
	endif
	if(highX < pnt2x(tempWave, 0)) //This does the replacing
		highX = pnt2x(tempWave, 0) 
	endif
	p1+=1
while(p1 < numWavesInList)
//Reset p1 again. Should have used For loops
p1=0
Do //Chops points off the end of the wave remembering that the deleted points are inclusive (thus I added one to the starting point).
	Wave tempWave = $(StringFromList(p1, inputWaveList) + "_faCorr")
	if(x2pnt(tempWave,lowX) > 0)
		DeletePoints 0,x2pnt(tempWave,highX) , tempWave
	endif
	Setscale/I x,highX,lowX, tempWave //This sets scale using the new starting X value and the old deltaX value (hence the /P). Has to be done here.
	p1+=1
while(p1 < numWavesInList)

Variable maxPntNum = 0 //Going to search for the maxium number of points in the list of waves just incase one of more waves has a drastically different deltaX
					     //Your deltaX should all be the same so you can probably ignore this bit.
//Reset p1 again out of habit. It's done automatically in the For loop.
p1 = 0
For(p1 = 0; p1< numWavesInList; p1 += 1)
	Wave tempWave = $(StringFromList(p1, inputWaveList)  + "_faCorr")
	If(numpnts(tempWave) > maxPntNum)
		maxPntNum = numpnts(tempWave)
	endif
endFor

//Print maxPntNum
//Reset p1 again out of habit. It's done automatically in the For loop.
p1 = 0
For(p1 = 0; p1< numWavesInList; p1 += 1) //Resample the waves to have the same number of points (ie fix any big differences in deltaX).
	Wave tempWave = $(StringFromList(p1, inputWaveList)  + "_faCorr")
	If(numpnts(tempWave) < maxPntNum)
		Interpolate2/N=(maxPntNum)/T=3/F=0 tempWave
	endif
EndFor
 
//Reset p1 again. Should have used For loops
p1 = 0
//Do //Undoing the Resampling I did at the begining and returning to as close to the original deltaX as possible. Usually get to within ~.0001 of the original value.
//	RatioFromNumber 1/newRate1
//	Wave tempWave = $(StringFromList(p1, inputWaveList) + "_faCorr")
//	Resample/UP=(V_numerator)/DOWN=(V_denominator)  tempWave
//	p1+=1
//while(p1 < numWavesInList)

end

Function WalkDataFoldersList1(path, waveListMatchStr, waveListOptionsStr)
	String path				// to walk all data folders, pass "root:" WITH the : char. To use the current data folder, use ":" or GetDataFolder(1)
	String waveListMatchStr		// "*" or "wave*" or "!*x", etc. See WaveList for details.
	String waveListOptionsStr	// usually "". see WaveList's optionsStr for details.
	SVAR outputWaveList = root:outputWaveList
	
	String savDF= GetDataFolder(1)
	SetDataFolder path
	
	//String listOfMatchingWaves= WaveList(waveListMatchStr,";",waveListOptionsStr)
	if(strlen(WaveList(waveListMatchStr,";",waveListOptionsStr)) && !stringMatch(GetDataFolder(0),"SensitivityFactors"))
		Variable p1 = 0
		for(p1 = 0; p1 < itemsInList(WaveList(waveListMatchStr,";",waveListOptionsStr)); p1+=1)
			Wave tempWave = $(StringFromList(p1,WaveList(waveListMatchStr,";",waveListOptionsStr)))
			outPutWaveList += (GetWavesDataFolder(tempWave,2) + ";")
		endfor
	endif
	
	Variable i, numDataFolders = CountObjects(":",4)
	for(i=0; i<numDataFolders; i+=1)
		String nextPath = GetIndexedObjName(":",4,i)
		WalkDataFoldersList1(nextPath, waveListMatchStr, waveListOptionsStr)
	endfor
	
	SetDataFolder savDF
End

Function WalkDataFoldersList(path, waveListMatchStr, waveListOptionsStr)
	String path				// to walk all data folders, pass "root:" WITH the : char. To use the current data folder, use ":" or GetDataFolder(1)
	String waveListMatchStr		// "*" or "wave*" or "!*x", etc. See WaveList for details.
	String waveListOptionsStr	// usually "". see WaveList's optionsStr for details.
	String/G root:outputWaveList = ""
	WalkDataFoldersList1(path, waveListMatchStr, waveListOptionsStr)
end

//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//
//																	 //
//																	 //
//  		          Two Pass First Derivative Guassian Convolution		        //
//																	 //
//																	 //
//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//

Function twoPassFDGMat(minWidth, maxWidth, numberPtsWidth, minEnergy, maxEnergy,  numberPtsSpec)
Variable minWidth, maxWidth, numberPtsWidth, minEnergy, maxEnergy, numberPtsSpec
Make/O/N=(numberPtsSpec) FDGauss
SetScale x,minEnergy,maxEnergy, FDGauss

Variable  tempWidth= abs(maxWidth-minWidth)/2, midPosition = abs(maxEnergy-minEnergy)/2

FDGauss =  -(x-midPosition)/(tempWidth^2)*exp(-(x-midPosition)^2/(2*tempWidth^2))

Make/O/N=(numberPtsSpec, numberPtsWidth) FDGaussMat
SetScale x,minEnergy,maxEnergy, FDGaussMat
SetScale/I y,minWidth,maxWidth, FDGaussMat
//Variable offSetCenter = -.5*(DimOffset(FDGaussMat, 1) + (numberPtsWidth-1)*DimDelta(FDGaussMat,1)) + pnt2x(FDGaussMat,0)
Variable offSetCenter = pnt2x(FDGaussMat,0) + ((ceil(numberPtsSpec/2))*deltaX(FDGaussMat))
//offSetCenter = pnt2x(FDGaussMat,0) - .5
FDGaussMat[][] =  -(x-offSetCenter)/(y^2)*exp(-(x-offSetCenter)^2/(2*y^2))
end

Function twoPassFDGconv(inputWave, minWidth, maxWidth, numberPtsWidth, percentHeight)
Wave inputWave
Variable minWidth, maxWidth, numberPtsWidth, percentHeight
Variable minEnergy = pnt2x(inputWave, 0), maxEnergy = pnt2x(inputWave, numpnts(inputWave)), numberPtsSpec = numpnts(inputWave)
String inputWaveName = NameOfWave(inputWave)

twoPassFDGMat(minWidth, maxWidth, numberPtsWidth, minEnergy, maxEnergy, numberPtsSpec)

Make/O/N=(numberPtsSpec,numberPtsWidth) $("convMat_" + inputWaveName)
Wave FDGaussMat, mat1 = $("convMat_" + inputWaveName)
SetScale x,minEnergy,maxEnergy, mat1
SetScale/I y,minWidth,maxWidth, mat1
Make/O/N=(numberPtsSpec) convTemp, temp
SetScale x,minEnergy,maxEnergy, convTemp

Variable p1 = 0, p2 = 0, p3 = 0
Do
	convTemp[] = FDGaussMat[p][p1]
	//MatrixOP/O temp = convolve(inputWave,convTemp,4)
	temp = inputWave
	Convolve/C convTemp,temp
	Convolve/C convTemp,temp
	mat1[][p1] = temp[p]
	p2 = 0
	Do
		If (mat1[p2][p1] < 0)
			mat1[p2][p1] = -mat1[p2][p1]
		else
			mat1[p2][p1] = 0
		endif
		p2 += 1
	while(p2 < dimsize(FDGaussMat,0))
	p1 += 1
while(p1 < dimsize(FDGaussMat,1))

localMaxSearch2D(mat1, percentHeight)
Wave localMaxWave
Make/D/O/N=(dimsize(localMaxWave,0),4) peakIntervals
Make/D/O/N=(dimsize(localMaxWave,0),2) peakSigmas

p1 = 0
Do
	p2 = 0
	Do
		p2 += 1
	while(mat1[localMaxWave[p1][2] + p2][localMaxWave[p1][3]] >= 10^-15 && (localMaxWave[p1][2] + p2) < (numpnts(inputWave)))
	p3 = 0
	Do
		p3 += 1
	while(mat1[localMaxWave[p1][2] - p3][localMaxWave[p1][3]] >= 10^-15 && (localMaxWave[p1][2] - p3) > 0)
	Variable sigmaL = 0, sigmaR = 0, xLeft, xRight
	xLeft = DimOffSet(mat1,0) + (localMaxWave[p1][2] - p3)*DimDelta(mat1,0)
	xRight = DimOffSet(mat1,0) + (localMaxWave[p1][2] + p2)*DimDelta(mat1,0)
	sigmaL = sqrt(abs((xLeft - localMaxWave[p1][0])^2 - 2*localMaxWave[p1][1]^2))
	sigmaR = sqrt(abs((xRight - localMaxWave[p1][0])^2 - 2*localMaxWave[p1][1]^2))
	peakSigmas[p1][0] = sigmaL
	peakSigmas[p1][1] = sigmaR
	peakIntervals[p1][0] = localMaxWave[p1][0] - 3*sigmaL
	peakIntervals[p1][1] = localMaxWave[p1][0] + 3*sigmaR
	peakIntervals[p1][2] = round((peakIntervals[p1][0] - DimOffSet(mat1,0))/DimDelta(mat1,0))
	peakIntervals[p1][3] = round((peakIntervals[p1][1] - DimOffSet(mat1,0))/DimDelta(mat1,0))
	If(peakIntervals[p1][2] < 0 || peakIntervals[p1][3] < 0)
	else
		Make/D/O/N=(abs(peakIntervals[p1][2] - peakIntervals[p1][3])) $("cutPeak_" + num2str(p1))
		Wave cutPeak = $("cutPeak_" + num2str(p1))
		cutPeak[] = inputWave[peakIntervals[p1][3] + p]
		SetScale x,pnt2x(inputWave, peakIntervals[p1][3]),pnt2x(inputWave, peakIntervals[p1][2]),cutPeak
	endif
	p1 += 1
while(p1 < dimsize(localMaxWave, 0))

Make/D/O/N=(numpnts(inputWave)) rWave
p1 = 0
p2 = 0

Do
	p2 = peakIntervals[p1][3]
	Do
		if(p2 < localMaxWave[p1][2])
			rWave[p2] = abs(peakIntervals[p1][3] - localMaxWave[p1][2])
		elseif(p2 >= localMaxWave[p1][2])
			rWave[p2] = abs(peakIntervals[p1][2] - localMaxWave[p1][2])
		endif
		p2 += 1
	while(p2 <= peakIntervals[p1][2])
	p1 += 1
while(p1 < dimsize(localMaxWave,0))

end

Function testDelta(inputWave, minWidth, maxWidth, numberPtsWidth)
Wave inputWave
Variable minWidth, maxWidth, numberPtsWidth
Variable deltaStart = 1, deltaEnd = numpnts(inputWave)
Variable numberPtsDelta = abs(deltaStart - deltaEnd), deltaDelta = (deltaEnd - deltaStart)/numberPtsDelta
Variable p1 = 0, deltaTemp = deltaStart
Make/O/N=(numpnts(inputWave),numberPtsWidth,numberPtsDelta) testDeltaMat

Do
	twoPassFDGconv(C1s, 0.03,1.02, 100, .3)
	Wave convMat_C1s
	testDeltaMat[][][p1] = convMat_C1s[p][q]
	deltaTemp += deltaDelta
	p1 += 1
while(p1 < numberPtsDelta)

end

Function markovChainSmth(inputWave, winSize)
Wave inputWave
Variable winSize
String inputWaveName = NameOfWave(inputWave)

Make/D/O/N=(numpnts(inputWave)) pPlus, pMinus, aWave ,$(inputWaveName + "_smth")
Wave smthWave = $(inputWaveName + "_smth")
SetScale x,pnt2x(inputWave,0),pnt2x(inputWave,numpnts(inputWave)),smthWave

smthWave = 0
pPlus = 0
pMinus = 0

Variable p1 = 0, p2 = 0, maxChan = WaveMax(inputWave)

Do
	p2 = 0
	Do
		pPlus[p1] += exp((inputWave[p1 + p2]/maxChan - inputWave[p1]/maxChan)/sqrt(inputWave[p1 + p2]/maxChan + inputWave[p1]/maxChan))
		pMinus[p1] += exp((inputWave[p1 - p2]/maxChan - inputWave[p1]/maxChan)/sqrt(inputWave[p1 - p2]/maxChan + inputWave[p1]/maxChan))
		p2 += 1
	While(p2 < winSize)
	p1 += 1
While(p1 < numpnts(inputWave))

aWave = 1/(pPlus + pMinus)
pPlus *= aWave
pMinus *= aWave

smthWave[0] = 1
p1 = 1
Do
	smthWave[p1] = (smthWave[p1 -1])*(pPlus[p1 - 1]/pMinus[p1])
	Variable sumWave = sum(smthWave)
	smthWave /= sumWave
	p1 += 1
While(p1 < numpnts(inputWave))
smthWave *= sum(inputWave)
smthWave += abs((wavemin(inputWave)-wavemin(smthWave)))
end

Function localMaxSearch2D(inputWave, percentHeight)
Wave inputWave
Variable percentHeight

Make/D/O/N=(0,4) localMaxWave

Variable p1 = 0, p2 = 0,p3 = 0 ,localMaxTemp = 0, localMaxflag = 0

Do
	p2 = 0
	Do
		mooreNeighborhood(inputWave,p1,p2)
		Wave mooreNeighbors
		p3 = 0
		Do
			If(mooreNeighbors[p3] > inputWave[p1][p2] || inputWave[p1][p2] <= 10^-15)
				localMaxflag = 0
				break
			elseif(p3 == numpnts(mooreNeighbors) - 1 && inputWave[p1][p2] >= percentHeight*Wavemax(inputWave))
				localMaxflag = 1
				InsertPoints/M=0 numpnts(localMaxWave),1,localMaxWave
				localMaxWave[numpnts(localMaxWave)][0] = DimOffset(inputWave,0) + p1*DimDelta(inputWave,0)
				localMaxWave[numpnts(localMaxWave)][1] = DimOffset(inputWave,1) + p2*DimDelta(inputWave,1)
				localMaxWave[numpnts(localMaxWave)][2] = p1
				localMaxWave[numpnts(localMaxWave)][3] = p2
			endif
			p3 += 1
		while(p3 < numpnts(mooreNeighbors))
		p2+=1
	while(p2 < dimsize(inputWave,1))
	p1+=1
while(p1 < dimsize(inputWave,0))

end

Function mooreNeighborhood(inputWave, xpnt, ypnt)
Wave inputWave
Variable xpnt, ypnt
Variable p1 = -1, p2 = -1, p3 = 0

Make/D/O/N=9 mooreNeighbors

Do
	p2 = -1
	Do
		mooreNeighbors[p3] = inputWave[xpnt + p1][ypnt + p2]
		p3 += 1
		p2 += 1
	while(p2 <= 1)
	p1 += 1
while(p1 <= 1)

DeletePoints 4,1,mooreNeighbors

end

//-------------------------------------------------------------------------------------------------------------//
// Adaptive Symmetric SNIP procedure									//
//-------------------------------------------------------------------------------------------------------------//

Function adaptSNIP(input_wave, rWave, order, smthWin)
Wave input_wave, rWave
Variable order, smthWin
String input_wave_str = NameOfWave(input_wave)
String nameHolder = input_wave_str + "_sub"
Make/O/N=(numpnts(input_wave)) workVect, baseline, LLSop_wave, $nameHolder
Wave w2 = $nameHolder 
Setscale x,pnt2x(input_wave,0),pnt2x(input_wave,numpnts(input_wave)),w2
Make/O/N=5 coefficients
Variable p1 = WaveMax(rWave),p2 = 0, p3 = 0, smthPnt = 0
//LLSop_wave = LLSop(input_wave_str)
LLSop_wave[] = input_wave[p]
Wave w1 = LLSop_wave
Do
	p2 = 0
	Do
	if(p1 <= rWave[p2])
		coefficients[0] = w1[p2]
		coefficients[1] = (w1[p2-p1]+w1[p2+p1])/2
		coefficients[2] = (-w1[p2-p1]+4*w1[p2-p1/2]+4*w1[p2+p1/2]-w1[p2+p1])/6
		coefficients[3] = (w1[p2-p1]-6*w1[p2-p1*2/3]+15*w1[p2-p1/3]+15*w1[p2+p1/3]-6*w1[p2+p1*2/3]+w1[p2+p1])/20
		coefficients[4] = (-w1[p2-p1]+8*w1[p2-p1*3/4]-28*w1[p2-p1/2]+56*w1[p2-p1/4]+56*w1[p2+p1/4]-28*w1[p2+p1/2]+8*w1[p2+p1*3/4]-w1[p2+p1])/70
		If(order == 2)
			workVect[p2] = min(coefficients[0],coefficients[1])
		elseif(order == 4)
			workVect[p2] = min(coefficients[0],max(coefficients[1],coefficients[2]))
		elseif(order == 6)
			workVect[p2] = min(coefficients[0],max(coefficients[1],max(coefficients[2],coefficients[3])))
		elseif(order == 8)
			workVect[p2] = min(coefficients[0],max(coefficients[1],max(coefficients[2],max(coefficients[3],coefficients[4]))))
		else
			break
		endif
	else
		p3 = -smthWin
		smthPnt = 0
		Do
			smthPnt += input_wave[p2 + p3]
			p3 +=1
		while(p3 <= smthWin)
		workVect[p2] = smthPnt/(2*smthWin+1)
	endif
		p2+=1
	while(p2 < numpnts(input_wave))
	w1 = workVect
	p1 -=1
while(p1 > 0)

//baseline = InverseLLSop("LLSop_wave")
baseline[] = w1[p]
WaveStats/Q input_wave
w2 = input_wave - baseline
end

Function XES_curve_calc(inputWave, slices)
Wave inputWave, slices
String inputWaveStr = NameOfWave(inputWave)

String maxiumList = ""
Variable p1 = 0
for(p1 = 0; p1 < dimsize(inputWave,1); p1+=1)
	If(slices[p1])
		Make/O/N=(dimsize(inputWave,0)) tempWave 
		tempWave[] = inputWave[p][p1]
		maxiumList +=  num2str(waveMaxLocation(tempWave)) + ";"
	endif
	KillWaves/Z tempWave
endfor

p1 = 0
Make/O/N=(itemsInList(maxiumList)) maxiumWave
for(p1 = 0; p1 < itemsInList(maxiumList); p1+=1)
	maxiumWave[p1] = str2num(stringfromList(p1, maxiumList, ";"))
endfor

SetScale/I x,3,28,maxiumWave
CurveFit/Q poly 3, maxiumWave/D
Wave fit_maxiumWave
Variable meanOfMax= fit_maxiumWave(dimsize(inputWave,1)/2 -1)

Make/O/N=(dimsize(inputWave,0),dimsize(inputWave,1)) $(inputWaveStr + "_curvCorr")
Wave outputWave = $(inputWaveStr + "_curvCorr"), W_coef
Variable p2 = 0
For(p1 = 0; p1 < dimsize(inputWave,0); p1 += 1)
	For(p2 = 0; p2 < dimsize(inputWave,1); p2 += 1)
		If(slices[p2])
			outputWave[p1][p2] = inputWave[p1+(W_coef[0]+ W_coef[1]*p2+W_coef[2]*p2^2)-meanOfMax][p2]
		else
			outputWave[p1][p2] = 0
		endIf
	endFor
endFor
Make/O/N=(dimsize(inputWave,0)) $(inputWaveStr + "_1D")
Wave outputWave1d = $(inputWaveStr + "_1D")
outputWave1d = 0
Variable sumCol = 0
For(p1 = 0; p1 < numpnts(outputWave1d); p1 += 1)
	For(p2 = 0; p2 < dimsize(inputWave,1); p2 += 1)
		If(slices[p2])
			outputWave1d[p1] += outputWave[p1][p2]
		endif
	endFor
endFor
KillWaves/Z maxiumWave,fit_maxiumWave, W_coef, W_sigma, W_ParamConfidenceInterval
end

Function waveMaxLocation(inputWave)
Wave inputWave
Duplicate/O/FREE inputwave,tempLocWave
//SetScale/P x,dimoffset(inputwave,0),dimdelta(inputwave,0),tempwave
Variable p1 = 0
Smooth/B=5 35, tempLocWave
Variable maximum = WaveMax(tempLocWave)
for(p1 = 0; p1 < numpnts(tempLocWave); p1+=1)
	If(tempLocWave[p1] == maximum)
		return pnt2x(tempLocWave,p1)
	endif
endfor
//FindPeak/B=5/M=(maximum*.8)/Q tempwave
//Return V_PeakLoc
end

Function autoShift_Gold(inputWaveList)
String inputWaveList
Variable p1 = 0
For(p1 = 0; p1 < itemsinlist(inputWaveList); p1 +=1)
	Wave tempWave = $stringfromlist(p1,inputWaveList)
	//Variable pntLocation = waveMaxLocation(tempWave)
	//pnt2x(tempWave,pntLocation)
	FindPeak/B=1/M=(WaveMax(tempWave)*.75) tempWave
	Variable xLocation = V_PeakLoc
	Variable photonEnergy = xLocation + 84
	Setscale/I x,(photonEnergy -pnt2x(tempWave,0)),(photonEnergy - pnt2x(tempWave,numpnts(tempWave))), tempWave 
endFor
end

Function autoShift_GoldVB(inputWaveList)
String inputWaveList
Variable p1 = 0
For(p1 = 0; p1 < itemsinlist(inputWaveList); p1 +=1)
	Wave tempWave = $stringfromlist(p1,inputWaveList)
	//Variable pntLocation = waveMaxLocation(tempWave)
	//pnt2x(tempWave,pntLocation)
	EdgeStats tempWave
	Variable xLocation = V_EdgeLoc2
	Variable photonEnergy = V_EdgeLoc2
	Setscale/I x,(photonEnergy -pnt2x(tempWave,0)),(photonEnergy - pnt2x(tempWave,numpnts(tempWave))), tempWave 
endFor
end

Function autoShift_toMax(inputWaveList, maxLoc)
String inputWaveList
Variable maxLoc
Variable p1 = 0
For(p1 = 0; p1 < itemsinlist(inputWaveList); p1 +=1)
	Wave tempWave = $stringfromlist(p1,inputWaveList)
	//FindPeak/M=(WaveMax(tempWave)*.9)/Q tempWave
	Variable xLocation = waveMaxLocation(tempWave)
	Variable shift = xLocation - maxLoc
	Setscale/I x,(pnt2x(tempWave,0)-shift),(pnt2x(tempWave,numpnts(tempWave))-shift), tempWave 
endFor
end

Function subMin(inputWaveList)
String inputWaveList
Variable p1 = 0
For(p1 = 0; p1 < itemsinlist(inputWaveList); p1 +=1)
	Wave tempWave = $stringfromlist(p1,inputWaveList)
	tempWave -= WaveMin(tempWave)
endFor
end

Function defineAtomicInfo()
String savDF = getDataFolder(1)
NewDataFolder/S root:AtomicInfo
Make/O  H_ ={1,1,1.0070,1}
Make/O  He_ ={2,3,4.0020,2}
Make/O  Li_ ={3,6,6.9410,1}
Make/O  Be_ ={4,9,9.0120,2}
Make/O  B_ ={5,10,10.8110,3}
Make/O  C_ ={6,12,12.0107,4}
Make/O  N_ ={7,14,14.0067,5}
Make/O  O_ ={8,16,15.9994,6}
Make/O  F_ ={9,19,18.9980,7}
Make/O  Ne_ ={10,20,20.1797,8}
Make/O  Na_ ={11,23,22.9890,1}
Make/O  Mg_ ={12,24,24.3050,2}
Make/O  Al_ ={13,27,26.9810,3}
Make/O  Si_ ={14,28,28.0855,4}
Make/O  P_ ={15,31,30.9730,5}
Make/O  S_ ={16,32,32.0650,6}
Make/O  Cl_ ={17,35,35.4530,7}
Make/O  Ar_ ={18,36,39.9480,8}
Make/O  K_ ={19,39,39.0983,1}
Make/O  Ca_ ={20,40,40.0780,2}
Make/O  Sc_ ={21,45,44.9550,3}
Make/O  Ti_ ={22,46,47.8670,4}
Make/O  V_ ={23,50,50.9415,5}
Make/O  Cr_ ={24,50,51.9961,6}
Make/O  Mn_ ={25,55,54.9380,7}
Make/O  Fe_ ={26,54,55.8450,8}
Make/O  Co_ ={27,59,58.9330,9}
Make/O  Ni_ ={28,58,58.6934,10}
Make/O  Cu_ ={29,63,63.5460,11}
Make/O  Zn_ ={30,64,65.3800,12}
Make/O  Ga_ ={31,69,69.7230,3}
Make/O  Ge_ ={32,70,72.6400,4}
Make/O  As_ ={33,75,74.9210,5}
Make/O  Se_ ={34,74,78.9600,6}
Make/O  Br_ ={35,79,79.9040,7}
Make/O  Kr_ ={36,78,83.7980,8}
Make/O  Rb_ ={37,85,85.4678,1}
Make/O  Sr_ ={38,84,87.6200,2}
Make/O  Y_ ={39,89,88.9050,3}
Make/O  Zr_ ={40,90,91.2240,4}
Make/O  Nb_ ={41,93,92.9060,5}
Make/O  Mo_ ={42,92,95.9600,6}
Make/O  Tc_ ={43,97,96.9060,7}
Make/O  Ru_ ={44,96,101.0700,8}
Make/O  Rh_ ={45,103,102.9050,9}
Make/O  Pd_ ={46,102,106.4200,10}
Make/O  Ag_ ={47,107,107.8682,11}
Make/O  Cd_ ={48,106,112.4110,12}
Make/O  In_ ={49,113,114.8180,3}
Make/O  Sn_ ={50,112,118.7100,4}
Make/O  Sb_ ={51,121,121.7600,5}
Make/O  Te_ ={52,120,127.6000,6}
Make/O  I_ ={53,127,126.9040,7}
Make/O  Xe_ ={54,124,131.2930,8}
Make/O  Cs_ ={55,133,132.9050,1}
Make/O  Ba_ ={56,130,137.3270,2}
Make/O  La_ ={57,138,138.9050,3}
Make/O  Ce_ ={58,136,140.1160,9}
Make/O  Pr_ ={59,141,140.9070,9}
Make/O  Nd_ ={60,142,144.2420,9}
Make/O  Pm_ ={61,145,144.9120,9}
Make/O  Sm_ ={62,144,150.3600,9}
Make/O  Eu_ ={63,151,151.9640,8}
Make/O  Gd_ ={64,152,157.2500,9}
Make/O  Tb_ ={65,159,158.9250,9}
Make/O  Dy_ ={66,156,162.5000,9}
Make/O  Ho_ ={67,165,164.9300,9}
Make/O  Er_ ={68,162,167.2590,9}
Make/O  Tm_ ={69,169,168.9340,9}
Make/O  Yb_ ={70,168,173.0540,8}
Make/O  Lu_ ={71,175,174.9668,9}
Make/O  Hf_ ={72,174,178.4900,4}
Make/O  Ta_ ={73,180,180.9470,5}
Make/O  W_ ={74,180,183.8400,6}
Make/O  Re_ ={75,185,186.2070,7}
Make/O  Os_ ={76,184,190.2300,8}
Make/O  Ir_ ={77,191,192.2170,9}
Make/O  Pt_ ={78,190,195.0840,10}
Make/O  Au_ ={79,197,196.9660,11}
Make/O  Hg_ ={80,196,200.5900,12}
Make/O  Tl_ ={81,203,204.3833,3}
Make/O  Pb_ ={82,204,207.2000,4}
Make/O  Bi_ ={83,209,208.9800,5}
Make/O  Po_ ={84,209,208.9820,6}
Make/O  At_ ={85,210,209.9870,7}
Make/O  Rn_ ={86,211,210.9900,8}
Make/O  Fr_ ={87,223,223.0190,1}
Make/O  Ra_ ={88,223,223.0180,2}
Make/O  Ac_ ={89,227,227.0270,3}
Make/O  Th_ ={90,230,230.0330,9}
Make/O  Pa_ ={91,231,231.0350,9}
Make/O  U_ ={92,233,233.0390,9}
Make/O  Np_ ={93,236,236.0460,9}
Make/O  Pu_ ={94,238,238.0490,9}
Make/O  Am_ ={95,241,241.0560,8}
Make/O  Cm_ ={96,243,243.0610,9}
Make/O  Bk_ ={97,247,247.0700,9}
Make/O  Cf_ ={98,249,249.0740,9}
Make/O  Es_ ={99,252,252.0820,9}
Make/O  Fm_ ={100,257,257.0950,9}
Make/O  Md_ ={101,258,258.0980,9}
Make/O  No_ ={102,259,259.1010,8}
Make/O  Lr_ ={103,262,262.1090,9}
Make/O  Rf_ ={104,265,265.1160,4}
Make/O  Db_ ={105,268,268.1250,5}
Make/O  Sg_ ={106,271,271.1330,6}
Make/O  Bh_ ={107,272,272.1380,7}
Make/O  Hs_ ={108,270,270.1340,8}
Make/O  Mt_ ={109,276,276.1510,9}
Make/O  Ds_ ={110,281,281.1620,10}
Make/O  Rg_ ={111,280,280.1640,11}
Make/O  Cn_ ={112,285,285.1740,12}
Make/O  Uut_ ={113,284,284.1780,3}
Make/O  Fl_ ={114,289,289.1870,4}
Make/O  Uup_ ={115,288,288.1920,5}
Make/O  Lv_ ={116,293,293.0000,6}
Make/O  Uus_ ={117,292,292.2070,7}
Make/O  Uuo_ ={118,294,294.0000,8}
SetDataFolder savDF
end

Function regionSelection(panelTop, panelBottom, panelLeft, panelRight, listPosX, listPosY, listSizeX, listSizeY)
Variable panelTop, panelBottom, panelLeft, panelRight, listPosX, listPosY, listSizeX, listSizeY
String/G root:Packages:regionSelectionReturnList = ""
String regionList, savDF = getDataFolder(1)

if(!DataFolderExists("root:SensitivityFactors"))
	defineSensitivityFactors()
endif

SetDataFolder root:SensitivityFactors
regionList = WaveList("!*_*",";","")
SetDataFolder savDF


Make/T/O/N=(itemsinlist(regionList)) root:Packages:listofRegions
Make/O/N=(itemsinList(regionList)) root:Packages:selectRegions

if(!WaveExists(root:Packages:listofSelectedWaves))
	Make/T/O/N=0 root:Packages:listOfSelectedWaves
	Make/O/N=0  root:Packages:deselectRegions
endif

Wave/T listofRegions = root:Packages:listofRegions, listOfSelectedWaves = root:Packages:listOfSelectedWaves
Wave selectRegions = root:Packages:selectRegions, deselectRegions = root:Packages:deselectRegions

Variable p1
For(p1 = 0; p1 < itemsinlist(regionList); p1 += 1)
	listofRegions[p1] = stringfromlist(p1, regionList)
endFor

DoWindow/K Region_Browser

NewPanel /W=(panelLeft,panelTop,panelRight,panelBottom)/N=Region_Browser as "Select Peak Regions"
ListBox regionSelectionBox,mode=9,pos={listPosX,listPosY},size={listSizeX, listSizeY},listWave= listofRegions,selWave=selectRegions
ListBox selectedRegionsBox,mode=9,pos={listPosX+listSizeX,listPosY},size={listSizeX, listSizeY},listWave= listofSelectedWaves,selWave=deselectRegions  
Button regionSelectionButton title="Select", pos={listPosX+10,listPosY+listSizeY+5}, size={80,20}, proc=regionSelectionButton_Proc
Button regionRemoveButton title="Remove",pos={listPosX+90,listPosY+listSizeY+5}, size={80,20}, proc=regionRemoveButton_Proc
Button regionSelectionDone title="Done",pos={listPosX+250,listPosY+listSizeY+5}, size={80,20}, proc=regionSelectCancel_Proc
Button regionSelectionCancel title="Cancel",pos={listPosX+170,listPosY+listSizeY+5}, size={80,20}, proc=regionSelectCancel_Proc
SetWindow Region_Browser, hook(regionSelectHookFunction_Cancel) = regionSelectHookFunction_Cancel
ModifyPanel fixedSize=1
end

//-------------------------------------------------------------------------------------------------------------------------------------------------------------------

Function regionSelectionButton_Proc(ctrlName) : ButtonControl
	String ctrlName
	SVAR regionSelectionReturnList = root:Packages:regionSelectionReturnList
	Wave/T listofRegions = root:Packages:listofRegions
	Wave selectRegions = root:Packages:selectRegions, inputAtomicMasses = root:Packages:inputAtomicMasses
	Wave/T listOfSelectedWaves = root:Packages:listOfSelectedWaves
	
	regionSelectionReturnList = ""
	Variable p1
	For(p1 = 0; p1 < numpnts(root:Packages:listOfSelectedWaves); p1 += 1)
		if(!stringMatch(listOfSelectedWaves[p1],""))
			regionSelectionReturnList = addListItem(listOfSelectedWaves[p1],regionSelectionReturnList,";",itemsinlist(regionSelectionReturnList))
		endif
	endFor
	
	For(p1 = 0; p1 < numpnts(listOfRegions); p1 += 1)
		if(selectRegions[p1] && strsearch(regionSelectionReturnList,listOfRegions[p1], 0) == -1)
			regionSelectionReturnList = addListItem(listOfRegions[p1],regionSelectionReturnList,";",itemsinList(regionSelectionReturnList))
		endif
	endfor
	Variable oldNumPnts = numpnts(root:Packages:listOfSelectedWaves)
	Redimension/N=(itemsInList(regionSelectionReturnList)) root:Packages:listOfSelectedWaves
	Wave/T listOfSelectedWaves = root:Packages:listOfSelectedWaves
	Redimension/N=(numpnts(root:Packages:listOfSelectedWaves)) root:Packages:deselectRegions, root:Packages:inputAtomicMasses
	Wave deselectRegions = root:Packages:deselectRegions
	deselectRegions = 0
	String tempListElements = ""
	For(p1 = 0; p1 < numpnts(root:Packages:listOfSelectedWaves); p1 += 1)
		listOfSelectedWaves[p1] = stringfromlist(p1, regionSelectionReturnList)
		String element 
		sscanf listOfSelectedWaves[p1],"%[A-Za-z]",element
		if(!stringMatch(element,""))
			if(strsearch(tempListElements,element,0) == -1)
				tempListElements = addlistitem(element,tempListElements,";",itemsinlist(tempListElements))
			else
				DoAlert/T="You Messed Up!" 0,"You can't have more than one region of any one element."
				Return 0
			endif
			Wave tempWave = $("root:AtomicInfo:" + element + "_")
			inputAtomicMasses[p1] = tempWave[2]
		endif
	endfor
	selectRegions = 0
End

Function regionRemoveButton_Proc(ctrlName) : ButtonControl
	String ctrlName
	SVAR regionSelectionReturnList = root:Packages:regionSelectionReturnList
	Wave/T listOfSelectedWaves = root:Packages:listOfSelectedWaves
	Wave deselectRegions = root:Packages:deselectRegions
	
	Variable p1
	For(p1 = 0; p1 < numpnts(deselectRegions); p1 += 1)
		if(deselectRegions[p1])
			regionSelectionReturnList = RemoveListItem(p1, regionSelectionReturnList)
			DeletePoints p1,1,listOfSelectedWaves,deselectRegions
			p1 -= 1
		endif
	endFor
	Make/O/N=(numpnts(listOfSelectedWaves))  root:Packages:deselectRegions = 0
	
End

//-------------------------------------------------------------------------------------------------------------------------------------------------------------------

Function regionSelectCancel_Proc(s) :ButtonControl
STRUCT WMButtonAction &s
Variable/G root:Packages:V_selectionResult
SVAR regionSelectionReturnList = root:Packages:regionSelectionReturnList
NVAR V_selectionResult = root:Packages:V_selectionResult
Wave/T listOfSelectedWaves = root:Packages:listOfSelectedWaves
//WS_ClearSelection("Wave_Browser", "waveSelectionBox")
if(s.eventCode == 2)
	if(stringMatch(s.ctrlName,"regionSelectionDone"))
		if(numpnts(listOfSelecetedWaves) > 15)
			V_selectionResult = 1
			DoAlert/T="You Messed Up!" 0,"There is no way you need that many regions. Try again."
			Return 0
		elseif(numpnts(listOfSelectedWaves) < 2)
			V_selectionResult = 1
			DoAlert/T="You Messed Up!" 0,"There has to be C1s and at least one other region. Try again."
			Return 0
		else
			Variable p1
			regionSelectionReturnList = ""
			For(p1 = 0; p1 < numpnts(listOfSelectedWaves); p1 += 1)
				regionSelectionReturnList = addListitem(listOfSelectedWaves[p1],regionSelectionReturnList,";",itemsInList(regionSelectionReturnList))
			endFor
			V_selectionResult = 0
		endif
	elseif(stringMatch(s.ctrlName,"regionSelectionCancel"))
		V_selectionResult = 1
		regionSelectionReturnList = ""
		Make/T/O/N=0 root:Packages:listofSelectedWaves
		Make/O/N=0 root:Packages:deselectRegions
		Make/T/O/N=0 root:Packages:regionArea_text, root:Packages:formula_text
		Make/O/N=0 root:Packages:areaHolderWave
		Make/O/N=0 root:Packages:formulaHolderWave, root:Packages:inputAtomicMasses
	endif
endif
genericCancel_Proc(s)

end

//-------------------------------------------------------------------------------------------------------------------------------------------------------------------

Function regionSelectHookFunction_Cancel(s)
STRUCT WMWinHookStruct &s
NVAR V_selectionResult = root:Packages:V_selectionResult
SVAR regionSelectionReturnList = root:Packages:regionSelectionReturnList
//WS_ClearSelection("Wave_Browser", "waveSelectionBox")
if(s.eventCode == 2)
	V_selectionResult = 1
	regionSelectionReturnList = ""
	Make/T/O/N=0 root:Packages:listofSelectedWaves
	Make/O/N=0 root:Packages:deselectRegions
	Make/T/O/N=0 root:Packages:regionArea_text
	Make/O/N=0 root:Packages:areaHolderWave
	KillWaves/Z root:Packages:listofRegions,root:Packages:selectRegions
endif
windowHookFunction_Cancel(s)

end

Function areaHookFunction_Cancel(s)
STRUCT WMWinHookStruct &s
Wave areaHolderWave = root:Packages:areaHolderWave
Make/T/O/N=(numpnts(areaHolderWave)) root:Packages:regionArea_text
Wave/T regionArea_text = root:Packages:regionArea_text
if(s.eventCode == 2)
	areaHolderWave = 0
		Variable p1
	For(p1 =0;p1 < numpnts(areaHolderWave);p1 +=1)
		regionArea_text[p1] = "0"
	endFor
endif
windowHookFunction_Cancel(s)

end

//-------------------------------------------------------------------------------------------------------------------------------------------------------------------

Function enterPeakAreas()
	SVAR regionSelectionReturnList =root:Packages:regionSelectionReturnList
	String dimensions = nearestDivisors(itemsinlist(regionSelectionReturnList))
	Variable largest, smallest
	if(str2num(stringfromlist(0,dimensions)) > str2num(stringfromlist(1,dimensions)))
		largest = str2num(stringfromlist(0,dimensions))
		smallest = str2num(stringfromlist(1,dimensions))
	else
		largest = str2num(stringfromlist(1,dimensions))
		smallest = str2num(stringfromlist(0,dimensions))
	endif
	Make/O/N=(itemsinlist(regionSelectionReturnList)) root:Packages:areaHolderWave
	Wave areaHolderWave = root:Packages:areaHolderWave
	NewPanel /W=(100,100,100+smallest*210+30,100+largest*20+50)/N=Peak_Areas as "Enter Peak Areas"
	Variable p1,p2,p3=0
	For(p1 = 0; p1 < largest; p1 += 1)
		For(p2 = 0; p2 < smallest; p2 += 1)
			SetVariable $("enterArea_" + num2str(p3)),win=Peak_Areas,title=stringfromlist(p3,regionSelectionReturnList),disable=0,limits={0,inf,0}
			SetVariable $("enterArea_" + num2str(p3)),win=Peak_Areas,size={200,20},pos={210*p2+15,20*p1+10},bodyWidth=150,value=areaHolderWave[p3]
			p3 += 1
		EndFor
	Endfor
	Button areaDoneButton,win=Peak_Areas,title="Done", pos={62,largest*20+15}, size={80,20},proc=areaEnterSelect_Proc
	Button areaCancelButton,win=Peak_Areas,title="Cancel",pos={142,largest*20+15}, size={80,20},proc=areaEnterCancel_Proc
	SetWindow Peak_Areas, hook(areaHookFunction_Cancel) = areaHookFunction_Cancel
	ModifyPanel fixedSize=1
end

Function/S divisors(inputNumber)
Variable inputNumber
PrimeFactors/Q inputNumber
Wave W_PrimeFactors
String factorList ="", outputList = ""

Variable p1, p2
For(p1 = 0; p1 < numpnts(W_PrimeFactors); p1 += 1)
	factorList += num2Str(W_PrimeFactors[p1]) + ";"
endFor
For(p1 = 0; p1 < itemsInlist(factorList); p1 += 1)
	For(p2 = p1; p2 < itemsInList(factorList); p2 +=1)
	if(mod(inputNumber,Str2num(stringfromlist(p1,factorList))*Str2num(stringfromlist(p2,factorList))) == 0 && Str2num(stringfromlist(p1,factorList))*Str2num(stringfromlist(p2,factorList)) < inputNumber && strsearch(factorList,num2Str(Str2num(stringfromlist(p1,factorList))*Str2num(stringfromlist(p2,factorList))),0) == -1)
		factorList += num2Str(Str2num(stringfromlist(p1,factorList))*Str2num(stringfromlist(p2,factorList)))+ ";"
	endif
	endFor
endFor

For(p1 = 0; p1 < itemsInList(factorList); p1 += 1)
	if(strSearch(outputList,stringFromlist(p1,factorList),0) == -1)
		outputList += stringFromlist(p1,factorList) +";"
	endif
endFor
Return Sortlist(outputList + "1;" + num2str(inputNumber) + ";", ";", 2)
end

Function/S nearestDivisors(inputNumber)
Variable inputNumber
Variable p1, tempDistance = inf
String factorList = divisors(inputNumber), tempPair = ""
For(p1 = 0; p1 < itemsinlist(factorList); p1 += 1)
	if(tempDistance > abs(inputNumber/str2num(stringfromlist(p1,factorList)) - str2num(stringfromlist(p1,factorList))))
		tempDistance = abs(inputNumber/str2num(stringfromlist(p1,factorList)) - str2num(stringfromlist(p1,factorList)))
		tempPair = stringfromlist(p1,factorList) +";"+num2str(inputNumber/str2num(stringfromlist(p1,factorList)))
	endif
endFor
Return tempPair
end

Function regionSelectStart_Proc(s) :ButtonControl
STRUCT WMButtonAction &s
if(s.eventCode == 2)
	regionSelection(100, 400, 100, 460, 10,10, 170, 250)
endif

end

Function areaEnterStart_Proc(s) :ButtonControl
STRUCT WMButtonAction &s
SVAR regionSelectionReturnList = root:Packages:regionSelectionReturnList
Wave listOfSelectedWaves = root:Packages:listOfSelectedWaves
if(s.eventCode == 2)
	DoWindow/K/Z Peak_Areas 
	if(numpnts(listOfSelectedWaves))
		enterPeakAreas()
	else
		DoAlert/T="You Messed Up!" 0,"There are no regions selected."
	endif
endif
end

Function areaEnterSelect_Proc(s) :ButtonControl
STRUCT WMButtonAction &s
Wave areaHolderWave = root:Packages:areaHolderWave
Make/T/O/N=(numpnts(areaHolderWave)) root:Packages:regionArea_text
Wave/T regionArea_text = root:Packages:regionArea_text
if(s.eventCode == 2)
	Variable p1
	For(p1 =0;p1 < numpnts(areaHolderWave);p1 +=1)
		regionArea_text[p1] = num2str(areaHolderWave[p1])
	endFor
endif
genericCancel_Proc(s)
end

Function areaEnterCancel_Proc(s) :ButtonControl
STRUCT WMButtonAction &s
Wave areaHolderWave = root:Packages:areaHolderWave
Make/T/O/N=(numpnts(areaHolderWave)) root:Packages:regionArea_text
Wave/T regionArea_text = root:Packages:regionArea_text
if(s.eventCode == 2)
	areaHolderWave = 0
		Variable p1
	For(p1 =0;p1 < numpnts(areaHolderWave);p1 +=1)
		regionArea_text[p1] = "0"
	endFor
endif
genericCancel_Proc(s)
end

Function OverLayerTab(s) : TabControl
	STRUCT WMTabControlAction &s
	
	ListBox selectedRegionList,disable= (s.tab!=0)
	TitleBox title_regionSelectionList,disable= (s.tab!=0)
	Button selectRegionsButton,disable= (s.tab!=0)
	ListBox areaList,disable= (s.tab!=0)
	TitleBox title_area,disable= (s.tab!=0)
	Button areasButton,disable= (s.tab!=0)
	ListBox formulaList,disable= (s.tab!=0)
	TitleBox title_formula,disable= (s.tab!=0)
	Button formulaButton,disable= (s.tab!=0)
	SetVariable distVal,disable= (s.tab!=0)
	SetVariable bandGapVal,disable= (s.tab!=0)
	SetVariable atomicWeightVal,disable= (s.tab!=0)
	SetVariable densityVal,disable= (s.tab!=0)
	SetVariable thetaVal,disable= (s.tab!=0)
	SetVariable photonEVal,disable= (s.tab!=0)
	Button calcButton,disable= (s.tab!=0)
	Button calcCancelButton,disable= (s.tab!=0)
End

Function iterativeOverLayerGUI()
	DoWindow/K/Z iterativeOverlayer
	NewDataFolder/O root:Packages
	if(!exists("root:Packages:regionSelectionReturnList"))
		Make/T/O/N=0 root:Packages:listofSelectedWaves, root:Packages:regionArea_text, root:Packages:formula_text
		Make/O/N=0 root:Packages:areaHolderWave, root:Packages:formulaHolderWave, root:Packages:deselectRegions
	else
		//SVAR regionSelectionReturnList = root:Packages:regionSelectionReturnList
		//Make/T/O/N=(itemsinlist(regionSelectionReturnList)) root:Packages:listofSelectedWaves, root:Packages:regionArea_text, root:Packages:formula_text
		//Make/O/N=(itemsinlist(regionSelectionReturnList)) root:Packages:areaHolderWave, root:Packages:formulaHolderWave
	endif
	Make/O/N=6 root:Packages:inputParameters = {15, 0, 0, 0, 45, 1486.7}
	Wave inputParameters = root:Packages:InputParameters
	NewPanel /K=1/W=(100,223,452,638)/N=iterativeOverlayer as "Iterative Overlayer Correction"
	ListBox selectedRegionList,pos={22,43},size={100,165}
	ListBox selectedRegionList,listWave=root:Packages:listOfSelectedWaves
	TitleBox title_regionSelectionList,pos={29,26},size={84,13},title="Selected Regions"
	TitleBox title_regionSelectionList,frame=0
	Button selectRegionsButton,pos={22,213},size={100,25},proc=regionSelectStart_Proc,title="Select Regions"
	ListBox areaList,pos={144,43},size={100,165}
	ListBox areaList,listWave=root:Packages:regionArea_text
	TitleBox title_area,pos={159,26},size={64,13},title="Region Areas",frame=0
	Button areasButton,pos={144,213},size={100,25},proc=areaEnterStart_Proc,title="Peak Areas"
	ListBox formulaList,pos={266,43},size={65,165}
	ListBox formulaList,listWave=root:Packages:formula_text
	TitleBox title_formula,pos={280,26},size={64,13},title="Formula",frame=0
	Button formulaButton,pos={266,213},size={65,25},proc=formulaEnterStart_Proc,title="Formula"
	SetVariable distVal,pos={86,247},size={206,16},bodyWidth=160,title="Distance"
	SetVariable distVal,limits={0,inf,0},value=inputParameters[0],help={"Initial Overlayer thickness in nm. Guess high."}
	SetVariable bandGapVal,pos={85,267},size={207,16},bodyWidth=160,title="Bandgap"
	SetVariable bandGapVal,limits={0,inf,0},value=inputParameters[1],help={"Bandgap in eV."}
	SetVariable atomicWeightVal,pos={56,287},size={236,16},bodyWidth=160,title="Atomic Weight "
	SetVariable atomicWeightVal,limits={0,inf,0},value=inputParameters[2],help={"Atomic Weight of a formula unit in amu."}
	SetVariable densityVal,pos={93,307},size={199,16},bodyWidth=160,title="Density"
	SetVariable densityVal,limits={0,inf,0},value=inputParameters[3],help={"Denisty in g/cm^3."}
	SetVariable thetaVal,pos={100,327},size={192,16},bodyWidth=160,title="Theta"
	SetVariable thetaVal,limits={0,inf,0},value=inputParameters[4],help={"Take of angle measured from the surface in degrees."}
	SetVariable photonEVal,pos={58,347},size={234,16},bodyWidth=160,title="Photon Energy"
	SetVariable photonEVal,limits={0,inf,0},value=inputParameters[5],help={"Photon Energy in eV. Usually shouldn't change this."}
	Button calcButton,pos={25,378},size={65,26},title="Calculate",proc=doIterativeCalc_Proc
	Button calcCancelButton,pos={95,378},size={65,26},title="Cancel",proc=iterativeCalcCancel_Proc
	TabControl overlayerTabs,pos={0,0},size={352,415},proc=OverLayerTab
	TabControl overlayerTabs,labelBack=(32768,32768,32768),tabLabel(0)="Setup"
	TabControl overlayerTabs,tabLabel(1)="Misc. Calc.",value= 0
	ModifyPanel fixedSize=1
end

Function iterativeCalcCancel_Proc(s) :ButtonControl
STRUCT WMButtonAction &s
genericCancel_Proc(s)
end

Function doIterativeCalc_Proc(s) :ButtonControl
STRUCT WMButtonAction &s
Wave/T listofSelectedWaves = root:Packages:listofSelectedWaves, formula_text = root:Packages:formula_text
Wave deselectRegions = root:Packages:deselectRegions, areaHolderWave = root:Packages:areaHolderWave, inputParameters =  root:Packages:inputParameters, formulaHolderWave = root:Packages:formulaHolderWave
Wave inputAtomicMasses = root:Packages:inputAtomicMasses
Variable p1
if(s.eventCode == 2)
	Make/O/N=(numpnts(listofSelectedWaves)) root:Packages:inputAtomicMasses
	Wave inputAtomicMasses = root:Packages:inputAtomicMasses
	Variable/G root:Packages:numValenceElec,root:Packages:initialCarbonArea
	NVAR numValenceElec = root:Packages:numValenceElec, intialCarbonArea = root:Packages:intialCarbonArea
	numValenceElec = 0
	String/G root:Packages:inputElementList  = ""
	SVAR inputElementList = root:Packages:inputElementList //Really the regions.
	String tempListElements = ""
	For(p1 = 0; p1 < numpnts(listofSelectedWaves); p1 += 1)
		inputElementList = addListItem(listofSelectedWaves[p1],inputElementList,";",itemsinlist(inputElementList))
		String element
		sscanf listOfSelectedWaves[p1],"%[A-Za-z]",element
		if(!stringMatch(element,""))
			if(strsearch(tempListElements,element,0) == -1)
				tempListElements = addlistitem(element,tempListElements,";",itemsinlist(tempListElements))
			else
				DoAlert/T="You Messed Up!" 0,"You can't have more than one region of any one element."
				Return 0
			endif
			Wave tempWave = $("root:AtomicInfo:" + element + "_")
			inputAtomicMasses[p1] = tempWave[2]
			numValenceElec += formulaHolderWave[p1]*tempWave[3]
		endif
	endFor
	
	if(!numpnts(listofSelectedWaves))
		DoAlert/T="You Messed Up!" 0,"You forgot to select some regions."
		Return 0
	elseif(!numpnts(areaHolderWave))
		DoAlert/T="You Messed Up!" 0,"You forgot to enter the peak areas."
		Return 0
	elseif(!numpnts(formulaHolderWave))
		DoAlert/T="You Messed Up!" 0,"You forgot to enter the chemical formula."
		Return 0
	elseif(strsearch(inputElementList,"C1s",0) == -1)
		DoAlert/T="You Messed Up!" 0,"Why are you trying to correct for a carbon overlayer if you don't have a C1s region?"
		Return 0
	endif
	For(p1 = 0; p1 < numpnts(inputParameters); p1 += 1)
		if(!inputParameters[p1])
			DoAlert/T="You Messed Up!" 0,"All input parameters must be nonzero."
			Return 0
		endif
	endFor
	For(p1 = 0; p1 < numpnts(areaHolderWave); p1 += 1)
		if(!areaHolderWave[p1])
			DoAlert/T="You Messed Up!" 0,"All areas must be nonzero."
			Return 0
		elseif(!formulaHolderWave[p1])
			if(strSearch(listofSelectedWaves[p1],"C1s",0) == -1)
				DoAlert/T="You Messed Up!" 0,"All non-carbon elements most have nonzero formula values."
				Return 0
			endif
		elseif(strSearch(listofSelectedWaves[p1],"C1s",0) != -1 && formulaHolderWave[p1])
				DoAlert/T="You Messed Up!" 0,"This procedure does not yet work for things containing carbon."
				Return 0
		endif
	endFor
	Variable C1sPsn = WhichListItem("C1s",inputElementList), C1sArea
	C1sArea = areaHolderWave[C1sPsn]
	inputElementList = RemoveListItem(C1sPsn,inputElementList)
	Make/O/N=(numpnts(listOfSelectedWaves) - 1) root:Packages:tempWaveAreas
	Wave tempWaveAreas = root:Packages:tempWaveAreas
	Variable p2 = 0
	For(p1 = 0; p1 < numpnts(tempWaveAreas)+1; p1 += 1)
		if(p1 != C1sPsn)
			tempWaveAreas[p2] = areaHolderWave[p1]
			p2 += 1
		endif
	endFor
	iterativeOverlayerCorrection(tempWaveAreas, tempWaveAreas , inputElementList, C1sArea, inputParameters[4], distance = inputParameters[0], density = inputParameters[3],  atomicWeight = inputParameters[2], BandGap = inputParameters[1], numValanceElectrons = numValenceElec, photonEnergy = inputParameters[5])
	Make/T/O/N=(numpnts(listOfSelectedWaves) - 1) regions
	Wave nValues
	Make/O/N=(numpnts(nValues)) percentComp
	percentComp = nValues*100
	For(p1 = 0; p1 < numpnts(regions); p1 += 1)
		regions[p1] = stringfromlist(p1,inputElementList)
	endFor
	DoWindow/K/Z Results
	Edit/K=1/N=Results regions,percentComp
endif
end

Function formulaEnterStart_Proc(s) :ButtonControl
STRUCT WMButtonAction &s
SVAR regionSelectionReturnList = root:Packages:regionSelectionReturnList
Wave listOfSelectedWaves = root:Packages:listOfSelectedWaves
if(s.eventCode == 2)
	DoWindow/K/Z Peak_Areas 
	if(numpnts(listOfSelectedWaves))
		enterFormula()
	else
		DoAlert/T="You Messed Up!" 0,"There are no regions selected."
	endif
endif
end


Function formulaEnterSelect_Proc(s) :ButtonControl
STRUCT WMButtonAction &s
Wave formulaHolderWave = root:Packages:formulaHolderWave
Make/T/O/N=(numpnts(formulaHolderWave)) root:Packages:formula_text
Wave/T formula_text = root:Packages:formula_text, listOfSelectedWaves = root:Packages:listOfSelectedWaves
Wave inputParameters = root:Packages:inputParameters, inputAtomicMasses = root:Packages:inputAtomicMasses

if(s.eventCode == 2)
	inputParameters[2] = 0
	Variable p1
	For(p1 =0;p1 < numpnts(formulaHolderWave);p1 += 1)
		formula_text[p1] = num2str(formulaHolderWave[p1])
	endFor
	For(p1 = 0; p1 < numpnts(formulaHolderWave); p1 += 1)
		if(!stringmatch(listOfSelectedWaves[p1],"C1s"))
			inputParameters[2] += formulaHolderWave[p1]*inputAtomicMasses[p1]
		endif
	endFor
endif
genericCancel_Proc(s)
end

Function formulaEnterCancel_Proc(s) :ButtonControl
STRUCT WMButtonAction &s
Wave formulaHolderWave = root:Packages:formulaHolderWave
Make/T/O/N=(numpnts(formulaHolderWave)) root:Packages:formula_text
Wave/T formula_text = root:Packages:formula_text
if(s.eventCode == 2)
	formulaHolderWave = 0
		Variable p1
	For(p1 =0;p1 < numpnts(formulaHolderWave);p1 +=1)
		formula_text[p1] = "0"
	endFor
endif
genericCancel_Proc(s)
end

Function enterFormula()
	SVAR regionSelectionReturnList =root:Packages:regionSelectionReturnList
	String dimensions = nearestDivisors(itemsinlist(regionSelectionReturnList))
	Variable largest, smallest
	if(str2num(stringfromlist(0,dimensions)) > str2num(stringfromlist(1,dimensions)))
		largest = str2num(stringfromlist(0,dimensions))
		smallest = str2num(stringfromlist(1,dimensions))
	else
		largest = str2num(stringfromlist(1,dimensions))
		smallest = str2num(stringfromlist(0,dimensions))
	endif
	Make/O/N=(itemsinlist(regionSelectionReturnList)) root:Packages:formulaHolderWave
	Wave formulaHolderWave = root:Packages:formulaHolderWave
	NewPanel /W=(100,100,100+smallest*210+30,100+largest*20+50)/N=Formula_ as "Chemical Formula"
	Variable p1,p2,p3=0
	For(p1 = 0; p1 < largest; p1 += 1)
		For(p2 = 0; p2 < smallest; p2 += 1)
			SetVariable $("enterArea_" + num2str(p3)),win=Formula_,title=stringfromlist(p3,regionSelectionReturnList),disable=0,limits={0,inf,0}
			SetVariable $("enterArea_" + num2str(p3)),win=Formula_,size={200,20},pos={210*p2+15,20*p1+10},bodyWidth=150,value=formulaHolderWave[p3]
			p3 += 1
		EndFor
	Endfor
	Button formulaDoneButton,win=Formula_,title="Done", pos={62,largest*20+15}, size={80,20},proc=formulaEnterSelect_Proc
	Button formulaCancelButton,win=Formula_,title="Cancel",pos={142,largest*20+15}, size={80,20},proc=formulaEnterCancel_Proc
	SetWindow Formula_, hook(areaHookFunction_Cancel) = areaHookFunction_Cancel
	ModifyPanel fixedSize=1
end

Function spinSplitAreaFactor(inputString)
String inputString
String element, orbital
sscanf inputString,"%[A-Za-z]%[A-Za-z0-9]",element,orbital
if(str2num(orbital[strlen(orbital)-1]))
	Variable jNum = str2num(orbital[strlen(orbital)-1])
	if(stringMatch(orbital[strlen(orbital)-2],"p"))
		switch(jNum)
			case 3:
				Return 3/2
				Break
			case 1:
				Return 3
				Break
		endswitch
	elseif(stringMatch(orbital[strlen(orbital)-2],"d"))
		switch(jNum)
			case 5:
				Return 5/3
				Break
			case 3:
				Return 5/2
				Break
		endswitch		
	elseif(stringMatch(orbital[strlen(orbital)-2],"f"))
		switch(jNum)
			case 7:
				Return 7/4
				Break
			case 5:
				Return 7/3
				Break
		endswitch		
	endif
else
	Return 1
endif

end

Function avePreEdge(inputWaveList)
String inputWaveList
Variable p1,p2,preEdgeAverage
For(p1 = 0; p1 < itemsInList(inputWaveList); p1 += 1)
	Wave tempWave = $(stringFromList(p1,inputWaveList)) 
	preEdgeAverage = 0
	For(p2 = numpnts(tempWave); p2 > numpnts(tempWave) - 20; p2 -= 1)
		preEdgeAverage += tempWave[p2]
	endFor
	preEdgeAverage /= 20
	tempWave -= preEdgeAverage
	Variable minVal = WaveMin(tempWave), maxVal = WaveMax(tempWave)
	tempWave = (tempWave - minVal)/(maxVal-minVal)
endFor
end

Function linearTheilSen(inputWave, xWave)
Wave inputWave, xWave
Variable numPairs = binomial(numpnts(inputWave),2)
Make /O/N=(numPairs) slopes
Make /O/N=(numPnts(inputWave)) intercepts
Variable p1,p2,count=0
For(p1 = 0; p1 < numpnts(inputWave); p1 += 1)
	For(p2 = p1+1; p2 < numpnts(inputWave); p2 += 1)
	slopes[count] = (inputWave[p1] - inputWave[p2])/(xWave[p1] - xWave[p2])
	count += 1
	endFor
endFor
Variable medianSlope = Median(slopes)
For(p1 = 0; p1 < numpnts(inputWave); p1 += 1)
	intercepts[p1] = inputWave[p1] - medianSlope*xWave[p1]
endFor
Variable medianIntercept = Median(intercepts)
Print/D medianSlope,medianIntercept
end

Function/D Median(w) // Returns median value of wave w
Wave w
Variable result
Duplicate/O w, temp // Make a clone of wave
Sort temp, temp // Sort clone
result = temp[numpnts(temp)/2]
KillWaves temp // Kill clone
return result
End

Function tricubic_spline_interp(inputVol,numInterpPnts)
Wave inputVol
Variable numInterpPnts

String inputVol_name = NameOfWave(inputVol)

Variable numNewPntsX, numNewPntsY, numNewPntsZ
numNewPntsX = dimsize(inputVol,0)*(numInterpPnts + 1) - numInterpPnts
numNewPntsY = dimsize(inputVol,1)*(numInterpPnts + 1) - numInterpPnts
numNewPntsZ = dimsize(inputVol,2)*(numInterpPnts + 1) - numInterpPnts
Make/O/N=(numNewPntsX,numNewPntsY,numNewPntsZ) $(inputVol_name + "_interp")
Wave outputVol = $(inputVol_name + "_interp")
outputVol = 0
Setscale/I x,DimOffset(inputVol,0),(dimsize(inputVol,0)-1)*DimDelta(inputVol,0)+DimOffset(inputVol,0),outputVol
Setscale/I y,DimOffset(inputVol,1),(dimsize(inputVol,1)-1)*DimDelta(inputVol,1)+DimOffset(inputVol,1),outputVol
Setscale/I z,DimOffset(inputVol,2),(dimsize(inputVol,2)-1)*DimDelta(inputVol,2)+DimOffset(inputVol,2),outputVol

Differentiate/DIM=0/METH=0  $(inputVol_name)/D=$(inputVol_name + "_xDiff")
Differentiate/DIM=1/METH=0  $(inputVol_name)/D=$(inputVol_name + "_yDiff")
Differentiate/DIM=2/METH=0  $(inputVol_name)/D=$(inputVol_name + "_zDiff")
Differentiate/DIM=1/METH=0  $(inputVol_name + "_xDiff")/D=$(inputVol_name + "_xyDiff")
Differentiate/DIM=2/METH=0  $(inputVol_name + "_xDiff")/D=$(inputVol_name + "_xzDiff")
Differentiate/DIM=2/METH=0  $(inputVol_name + "_yDiff")/D=$(inputVol_name + "_yzDiff")
Differentiate/DIM=2/METH=0  $(inputVol_name + "_xyDiff")/D=$(inputVol_name + "_xyzDiff")

Make/O/N=(64,64) bMat_inv 
    bMat_inv [][0] = { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][1] = { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][2] = {-3, 3, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][3] = { 2,-2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][4] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][5] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][6] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][7] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][8] = {-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][9] = { 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][10] = { 9,-9,-9, 9, 0, 0, 0, 0, 6, 3,-6,-3, 0, 0, 0, 0, 6,-6, 3,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][11] = {-6, 6, 6,-6, 0, 0, 0, 0,-3,-3, 3, 3, 0, 0, 0, 0,-4, 4,-2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-2,-1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][12] = { 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][13] = { 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][14] = {-6, 6, 6,-6, 0, 0, 0, 0,-4,-2, 4, 2, 0, 0, 0, 0,-3, 3,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][15] = { 4,-4,-4, 4, 0, 0, 0, 0, 2, 2,-2,-2, 0, 0, 0, 0, 2,-2, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][16] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][17] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][18] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][19] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][20] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][21] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][22] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0}
    bMat_inv [][23] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0}
    bMat_inv [][24] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][25] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0}
    bMat_inv [][26] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9,-9,-9, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 3,-6,-3, 0, 0, 0, 0, 6,-6, 3,-3, 0, 0, 0, 0, 4, 2, 2, 1, 0, 0, 0, 0}
    bMat_inv [][27] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3,-3, 3, 3, 0, 0, 0, 0,-4, 4,-2, 2, 0, 0, 0, 0,-2,-2,-1,-1, 0, 0, 0, 0}
    bMat_inv [][28] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][29] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0}
    bMat_inv [][30] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4,-2, 4, 2, 0, 0, 0, 0,-3, 3,-3, 3, 0, 0, 0, 0,-2,-1,-2,-1, 0, 0, 0, 0}
    bMat_inv [][31] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4,-4,-4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2,-2,-2, 0, 0, 0, 0, 2,-2, 2,-2, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0}
    bMat_inv [][32] = {-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][33] = { 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][34] = { 9,-9, 0, 0,-9, 9, 0, 0, 6, 3, 0, 0,-6,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6,-6, 0, 0, 3,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 2, 0, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][35] = {-6, 6, 0, 0, 6,-6, 0, 0,-3,-3, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 4, 0, 0,-2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-2, 0, 0,-1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][36] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][37] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0}
    bMat_inv [][38] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9,-9, 0, 0,-9, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 3, 0, 0,-6,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6,-6, 0, 0, 3,-3, 0, 0, 4, 2, 0, 0, 2, 1, 0, 0}
    bMat_inv [][39] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 0, 0, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3,-3, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 4, 0, 0,-2, 2, 0, 0,-2,-2, 0, 0,-1,-1, 0, 0}
    bMat_inv [][40] = { 9, 0,-9, 0,-9, 0, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 3, 0,-6, 0,-3, 0, 6, 0,-6, 0, 3, 0,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][41] = { 0, 0, 0, 0, 0, 0, 0, 0, 9, 0,-9, 0,-9, 0, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 3, 0,-6, 0,-3, 0, 6, 0,-6, 0, 3, 0,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 2, 0, 1, 0}
    bMat_inv [][42] = {-27,27,27,-27,27,-27,-27,27,-18,-9,18, 9,18, 9,-18,-9,-18,18,-9, 9,18,-18, 9,-9,-18,18,18,-18,-9, 9, 9,-9,-12,-6,-6,-3,12, 6, 6, 3,-12,-6,12, 6,-6,-3, 6, 3,-12,12,-6, 6,-6, 6,-3, 3,-8,-4,-4,-2,-4,-2,-2,-1}
    bMat_inv [][43] = {18,-18,-18,18,-18,18,18,-18, 9, 9,-9,-9,-9,-9, 9, 9,12,-12, 6,-6,-12,12,-6, 6,12,-12,-12,12, 6,-6,-6, 6, 6, 6, 3, 3,-6,-6,-3,-3, 6, 6,-6,-6, 3, 3,-3,-3, 8,-8, 4,-4, 4,-4, 2,-2, 4, 4, 2, 2, 2, 2, 1, 1}
    bMat_inv [][44] = {-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0,-3, 0, 3, 0, 3, 0,-4, 0, 4, 0,-2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-2, 0,-1, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][45] = { 0, 0, 0, 0, 0, 0, 0, 0,-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0,-3, 0, 3, 0, 3, 0,-4, 0, 4, 0,-2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-2, 0,-1, 0,-1, 0}
    bMat_inv [][46] = {18,-18,-18,18,-18,18,18,-18,12, 6,-12,-6,-12,-6,12, 6, 9,-9, 9,-9,-9, 9,-9, 9,12,-12,-12,12, 6,-6,-6, 6, 6, 3, 6, 3,-6,-3,-6,-3, 8, 4,-8,-4, 4, 2,-4,-2, 6,-6, 6,-6, 3,-3, 3,-3, 4, 2, 4, 2, 2, 1, 2, 1}
    bMat_inv [][47] = {-12,12,12,-12,12,-12,-12,12,-6,-6, 6, 6, 6, 6,-6,-6,-6, 6,-6, 6, 6,-6, 6,-6,-8, 8, 8,-8,-4, 4, 4,-4,-3,-3,-3,-3, 3, 3, 3, 3,-4,-4, 4, 4,-2,-2, 2, 2,-4, 4,-4, 4,-2, 2,-2, 2,-2,-2,-2,-2,-1,-1,-1,-1}
    bMat_inv [][48] = { 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][49] = { 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][50] = {-6, 6, 0, 0, 6,-6, 0, 0,-4,-2, 0, 0, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][51] = { 4,-4, 0, 0,-4, 4, 0, 0, 2, 2, 0, 0,-2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][52] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][53] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0}
    bMat_inv [][54] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 0, 0, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4,-2, 0, 0, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-3, 3, 0, 0,-2,-1, 0, 0,-2,-1, 0, 0}
    bMat_inv [][55] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4,-4, 0, 0,-4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0,-2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 2,-2, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0}
    bMat_inv [][56] = {-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 0,-2, 0, 4, 0, 2, 0,-3, 0, 3, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][57] = { 0, 0, 0, 0, 0, 0, 0, 0,-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 0,-2, 0, 4, 0, 2, 0,-3, 0, 3, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0,-2, 0,-1, 0}
    bMat_inv [][58] = {18,-18,-18,18,-18,18,18,-18,12, 6,-12,-6,-12,-6,12, 6,12,-12, 6,-6,-12,12,-6, 6, 9,-9,-9, 9, 9,-9,-9, 9, 8, 4, 4, 2,-8,-4,-4,-2, 6, 3,-6,-3, 6, 3,-6,-3, 6,-6, 3,-3, 6,-6, 3,-3, 4, 2, 2, 1, 4, 2, 2, 1}
    bMat_inv [][59] = {-12,12,12,-12,12,-12,-12,12,-6,-6, 6, 6, 6, 6,-6,-6,-8, 8,-4, 4, 8,-8, 4,-4,-6, 6, 6,-6,-6, 6, 6,-6,-4,-4,-2,-2, 4, 4, 2, 2,-3,-3, 3, 3,-3,-3, 3, 3,-4, 4,-2, 2,-4, 4,-2, 2,-2,-2,-1,-1,-2,-2,-1,-1}
    bMat_inv [][60] = { 4, 0,-4, 0,-4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 2, 0,-2, 0,-2, 0, 2, 0,-2, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    bMat_inv [][61] = { 0, 0, 0, 0, 0, 0, 0, 0, 4, 0,-4, 0,-4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 2, 0,-2, 0,-2, 0, 2, 0,-2, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0}
    bMat_inv [][62] = {-12,12,12,-12,12,-12,-12,12,-8,-4, 8, 4, 8, 4,-8,-4,-6, 6,-6, 6, 6,-6, 6,-6,-6, 6, 6,-6,-6, 6, 6,-6,-4,-2,-4,-2, 4, 2, 4, 2,-4,-2, 4, 2,-4,-2, 4, 2,-3, 3,-3, 3,-3, 3,-3, 3,-2,-1,-2,-1,-2,-1,-2,-1}
    bMat_inv [][63] = { 8,-8,-8, 8,-8, 8, 8,-8, 4, 4,-4,-4,-4,-4, 4, 4, 4,-4, 4,-4,-4, 4,-4, 4, 4,-4,-4, 4, 4,-4,-4, 4, 2, 2, 2, 2,-2,-2,-2,-2, 2, 2,-2,-2, 2, 2,-2,-2, 2,-2, 2,-2, 2,-2, 2,-2, 1, 1, 1, 1, 1, 1, 1, 1}
    
MatrixTranspose bMat_inv
Make/O/N=64 coefficients  = 0

Make/O/N=64 $(inputVol_name + "_constraints")
Wave constraints = $(inputVol_name + "_constraints")

Variable p0,p1,p2,c0,n0,n1,n2,m0,m1,m2,p0_new,p1_new,p2_new,interpValue,k0,k1,k2,x_out,y_out,z_out
String waveString

For(p0 = 0; p0 < dimsize(inputVol,0) -1; p0 += 1)
	For(p1 = 0; p1 < dimsize(inputVol,1) -1; p1 += 1)
		For(p2 = 0; p2 < dimsize(inputVol,2) -1; p2 += 1)
			For(c0 = 0; c0 < 64; c0 += 8)
				If(c0 == 0)
					waveString = inputVol_name
				Elseif(c0 == 8)
					waveString = inputVol_name + "_xDiff"
				Elseif(c0 == 16)
					waveString = inputVol_name + "_yDiff"
				Elseif(c0 == 24)
					waveString = inputVol_name + "_zDiff"
				Elseif(c0 == 32)
					waveString = inputVol_name + "_xyDiff"
				Elseif(c0 == 40)
					waveString = inputVol_name + "_xzDiff"
				Elseif(c0 == 48)
					waveString = inputVol_name + "_yzDiff"
				Elseif(c0 == 56)
					waveString = inputVol_name + "_xyzDiff"
				endIf
				Wave currentWave = $waveString
				constraints[c0+0] = currentWave[p0][p1][p2]
				constraints[c0+1] = currentWave[p0+1][p1][p2]
				constraints[c0+2] = currentWave[p0][p1+1][p2]
				constraints[c0+3] = currentWave[p0+1][p1+1][p2]
				constraints[c0+4] = currentWave[p0][p1][p2+1]
				constraints[c0+5] = currentWave[p0+1][p1][p2+1]
				constraints[c0+6] = currentWave[p0][p1+1][p2+1]
				constraints[c0+7] = currentWave[p0+1][p1+1][p2+1]
			endFor
			
			MatrixMultiply bMat_inv,constraints
			Wave M_product
			coefficients = M_product
			
			p0_new = p0*(numInterpPnts+1)
			p1_new = p1*(numInterpPnts+1)
			p2_new = p2*(numInterpPnts+1)
			
			For(m0 = 0; m0 <= numInterpPnts + 1; m0 += 1)
				For(m1 = 0; m1 <= numInterpPnts + 1 ; m1 += 1)
					For(m2 = 0; m2 <= numInterpPnts + 1; m2 += 1)
						x_out = m0/(numInterpPnts+1)
						y_out = m1/(numInterpPnts+1)
						z_out = m2/(numInterpPnts+1)
						For(n0 = 0; n0 < 4; n0 += 1)
							For(n1 = 0; n1 < 4; n1 += 1)
								For(n2 = 0; n2 < 4; n2 += 1)
									interpValue += coefficients[n0+4*n1+16*n2]*(x_out^(n0))*(y_out^(n1))*(z_out^(n2))
								endFor
							endFor
						endFor
						outputVol[p0_new+m0][p1_new+m1][p2_new+m2] = interpValue
						interpvalue = 0
					endFor
				endFor
			endFor
		endFor
	endFor
endFor	
end

Function fd_coeff(m)
Variable m
Variable j,k
Make/O/N=(2*m+1,2*m+1) fd_coeff_mat

For(j = 0; j < 5; j += 1)
	For(k = 0; k < 5; k += 1)
		fd_coeff_mat[k][j] = ((k+1) - m - 1)^((j+1)-1)/factorial((j+1) - 1)
	endFor
endFor

MatrixInverse fd_coeff_mat
end

Function fd_derivative(inputWave,m)
Wave inputWave
Variable m
Duplicate/O inputWave, $(NameOfWave(inputWave) + "_Diff")
Wave outWave = $(NameOfWave(inputWave) + "_Diff")
fd_coeff(m)
Wave fd_coeff_mat_inv = M_inverse
Variable p0,p1,result
For(p0 = 0; p0 < numpnts(inputWave); p0 += 1)
	result = 0
	For(p1 = -m; p1 <= m; p1 += 1)
		result += fd_coeff_mat_inv[1][p1 + m]*inputWave[p0+p1]
	endFor
	outWave[p0] = result
endFor
end

Function LCASI_1D_periodic(inputWave, degree)
Wave inputWave
Variable degree

Variable i,j,interp_i,x_value,left_window,right_window
Make/N=3 /FREE y_left,y_right,coeff_left, coeff_right
Make/N=(3,3) /FREE inv_left = {{0.5,-0.5,0.0},{-1.0,0.0,1.0},{0.5,0.5,0.0}},inv_right = {{0.5,-1.5,1.0},{-1.0,2.0,0.0},{0.5,-0.5,0.0}}

Variable numPnts_input = numPnts(inputWave),numPnts_interp = (numpnts(inputWave)-1)*degree+1

Make/FREE/N=(numPnts_interp) outWave
//Wave outWave = $(NameOfWave(inputWave) + "_interp")
setScale/P x,DimOffset(inputWave,0),DimDelta(inputWave,0)/degree,outWave

y_left[0] = inputWave[numPnts_input - 1]
y_left[1] = inputWave[0]
y_left[2] = inputWave[1]

MatrixOp/O coeff_left = inv_left x y_left
For(i = 0; i < numPnts_input; i += 1)
	interp_i = i*degree
	For(j = 0; j < 3; j += 1)
		if(i + j < numPnts_input)
			y_right[j] = inputWave[i+j]
		else
			y_right[j] = inputWave[i+j-numPnts_input]
		endif
	endFor
	MatrixOp/O coeff_right = inv_right x y_right
	outWave[interp_i] = inputWave[i]
	For (j = 1; j < degree; j += 1)
		x_value = j / degree
		right_window = x_value
		left_window = 1.0 - x_value
		if(interp_i+j < numPnts_interp)
			outWave[interp_i+j] = left_window*(x_value*(coeff_left[0]*x_value + coeff_left[1]) + coeff_left[2]) + right_window*(x_value*(coeff_right[0]*x_value + coeff_right[1]) + coeff_right[2])
		endif
	endFor
	For(j = 0; j < 3; j += 1)
		if(i + j < numPnts_input)
			y_left[j] = inputWave[i+j]
		else
			y_left[j] = inputWave[i+j-numPnts_input]
		endif
	endFor
	MatrixOp/O coeff_left = inv_left x y_left	
endFor
	return outWave
end

Function LCASI_2D_periodic(inputWave, outputWave, degree_x, degree_y)
Wave inputWave, outputWave
Variable degree_x, degree_y

Variable numPnts_input_x = dimSize(inputWave,0),numPnts_input_y = dimSize(inputWave,1)
Variable numPnts_interp_x = (dimSize(inputWave,0)-1)*degree_x+1,numPnts_interp_y = (dimSize(inputWave,1)-1)*degree_y+1

Make/FREE/N=(numPnts_interp_x,numPnts_interp_y) outWave_2D
setScale/P x,DimOffset(inputWave,0),DimDelta(inputWave,0)/degree_x,outWave_2D
setScale/P y,DimOffset(inputWave,1),DimDelta(inputWave,1)/degree_y,outWave_2D

Variable p0,p1
For(p0 = 0; p0 < numPnts_input_x; p0 += 1)
	MatrixOp/FREE inputWaveSection = row(inputWave,p0)^t
	Make/FREE/N=(numPnts_interp_y) outWaveSection = LCASI_1D_periodic(inputWaveSection, degree_y)
	outWave_2D[p0*degree_x][] = outWaveSection[q]
endFor
For(p1 = 0; p1 < numPnts_interp_y; p1 += 1)
	Make/FREE/N=(numPnts_input_x) inputWaveSection = outWave_2D[p*degree_x][p1]
	Make/FREE/N=(numPnts_interp_x) outWaveSection = LCASI_1D_periodic(inputWaveSection, degree_x)
	outWave_2D[][p1] = outWaveSection[p]		
endFor
outputWave[][] = outWave_2D[p][q]
end

Function LCASI_3D_periodic(inputWave, outputWave, degree_x, degree_y, degree_z)
Wave inputWave, outputWave
Variable degree_x,degree_y,degree_z

Variable numPnts_input_x = dimSize(inputWave,0),numPnts_input_y = dimSize(inputWave,1),numPnts_input_z = dimSize(inputWave,2)
Variable numPnts_interp_x = (dimSize(inputWave,0)-1)*degree_x+1,numPnts_interp_y = (dimSize(inputWave,1)-1)*degree_y+1,numPnts_interp_z = (dimSize(inputWave,2)-1)*degree_z+1

Make/FREE/N=(numPnts_interp_x,numPnts_interp_y,numPnts_interp_z) outWave_3D
setScale/P x,DimOffset(inputWave,0),DimDelta(inputWave,0)/degree_x,outWave_3D
setScale/P y,DimOffset(inputWave,1),DimDelta(inputWave,1)/degree_y,outWave_3D
setScale/P z,DimOffset(inputWave,2),DimDelta(inputWave,2)/degree_z,outWave_3D

Variable p0,p1
For(p0 = 0; p0 < numPnts_input_z; p0 += 1)
	Make/FREE/N=(numPnts_input_x,numPnts_input_y) inputWaveSection = inputWave[p][q][p0]
	Make/FREE/N=(numPnts_interp_x,numPnts_interp_y) outWaveSection
	LCASI_2D_periodic(inputWaveSection, outWaveSection, degree_x, degree_y)
	outWave_3D[][][p0*degree_z] = outWaveSection[p][q]
endFor
For(p1 = 0; p1 < numPnts_interp_x; p1 += 1)
	Make/FREE/N=(numPnts_input_y,numPnts_input_z) inputWaveSection = outWave_3D[p1][p*degree_y][q*degree_z]
	Make/FREE/N=(numPnts_interp_y,numPnts_interp_z) outWaveSection
	LCASI_2D_periodic(inputWaveSection, outWaveSection, degree_y, degree_z)
	outWave_3D[p1][][] = outWaveSection[q][r]
endFor

outputWave[][][] = outWave_3D[p][q][r]
end

Function Interp3D_LCASI(inputWave, degree_x, degree_y, degree_z)
Wave inputWave
Variable degree_x,degree_y,degree_z

Variable numPnts_input_x = dimSize(inputWave,0),numPnts_input_y = dimSize(inputWave,1),numPnts_input_z = dimSize(inputWave,2)
Variable numPnts_interp_x = (dimSize(inputWave,0)-1)*degree_x+1,numPnts_interp_y = (dimSize(inputWave,1)-1)*degree_y+1,numPnts_interp_z = (dimSize(inputWave,2)-1)*degree_z+1

Make/O/N=(numPnts_interp_x,numPnts_interp_y,numPnts_interp_z) $(nameOfWave(inputWave) + "_interp")
Wave outputWave_3D = $(nameOfWave(inputWave) + "_interp")
setScale/P x,DimOffset(inputWave,0),DimDelta(inputWave,0)/degree_x,outputWave_3D
setScale/P y,DimOffset(inputWave,1),DimDelta(inputWave,1)/degree_y,outputWave_3D
setScale/P z,DimOffset(inputWave,2),DimDelta(inputWave,2)/degree_z,outputWave_3D
LCASI_3D_periodic(inputWave, outputWave_3D, degree_x, degree_y, degree_z)
end

Function Interp2D_LCASI(inputWave, degree_x, degree_y)
Wave inputWave
Variable degree_x,degree_y

Variable numPnts_input_x = dimSize(inputWave,0),numPnts_input_y = dimSize(inputWave,1)
Variable numPnts_interp_x = (dimSize(inputWave,0)-1)*degree_x+1,numPnts_interp_y = (dimSize(inputWave,1)-1)*degree_y+1

Make/O/N=(numPnts_interp_x,numPnts_interp_y) $(nameOfWave(inputWave) + "_interp")
Wave outputWave_2D = $(nameOfWave(inputWave) + "_interp")
setScale/P x,DimOffset(inputWave,0),DimDelta(inputWave,0)/degree_x,outputWave_2D
setScale/P y,DimOffset(inputWave,1),DimDelta(inputWave,1)/degree_y,outputWave_2D
LCASI_2D_periodic(inputWave, outputWave_2D, degree_x, degree_y)
end

Function Interp1D_LCASI(inputWave, degree_x)
Wave inputWave
Variable degree_x

Variable numPnts_input_x = dimSize(inputWave,0)
Variable numPnts_interp_x = (dimSize(inputWave,0)-1)*degree_x+1

Make/O/N=(numPnts_interp_x) $(nameOfWave(inputWave) + "_interp")
Wave outputWave_1D = $(nameOfWave(inputWave) + "_interp")
setScale/P x,DimOffset(inputWave,0),DimDelta(inputWave,0)/degree_x,outputWave_1D
outputWave_1D = LCASI_1D_periodic(inputWave, degree_x)
end

Function sumX24a(inputWave)
Wave inputWave
String inputWaveStr = NameOfWave(inputWave)
String waveNote = note(inputWave)

Make/O/N=(dimsize(inputWave,0)) $(inputWaveStr + "_1D") = 0
Wave outWave = $(inputWaveStr + "_1D")
Note outWave, waveNote

Setscale/P x,dimoffset(inputWave,0),dimdelta(inputWave,0), outWave
Variable p1 = 0,p2 = 0
Do
	Do
		outWave[] += inputWave[p][p2][p1]
		p2 += 1
	While(p2 < dimsize(inputWave,1))
	p1 += 1
While(p1 < dimsize(inputWave,2))
end

Function beamDriftCor(inputMat)
Wave inputMat
If(dimsize(inputMat,1))
	Make/O/N=(dimsize(inputMat,0),dimsize(inputMat,1))/FREE tempMat = inputMat
	Variable p1 = 0
	String listOfWaves = ""
	DFREF dfrSave = GetDataFolderDFR()
	DFREF dfrFree =  NewFreeDataFolder()
	SetDataFolder dfrFree

	For(p1 = 0; p1 < dimsize(tempMat,1); p1 += 1)
		Make/O/N=(dimsize(inputMat,0)) $("tempWave_" + num2str(p1)) = inputMat[p][p1]
		SetScale/P x,dimoffset(inputMat,0),dimdelta(inputMat,0),$("tempWave_" + num2str(p1))
		Interpolate2/Y=$("tempWave_" + num2str(p1) + "_interp")/T=1/N=(dimSize(inputMat,0)*10) $("tempWave_" + num2str(p1))
		listOfWaves = addListItem("tempWave_" + num2str(p1) + "_interp",listOfWaves,";",inf)
	Endfor
	autoShift_toMax(listOfWaves,waveMaxLocation(tempWave_0_interp))
	WaveShave(listOfWaves)
	Duplicate/O inputMat, outputMat
	Redimension/N=(dimsize(tempWave_0_interp_faCorr,0),dimsize(outputMat,1)) outputMat
	SetScale/P x,dimoffset(tempWave_0_interp_faCorr,0),dimdelta(tempWave_0_interp_faCorr,0),outputMat

	For(p1 = 0; p1 < dimsize(tempMat,1); p1 += 1)
		wave holderWave = $(stringFromList(p1, listOfWaves) + "_faCorr")
		outputMat[][p1] = holderWave[p]
	EndFor
	SetDataFolder dfrSave
	//MoveDataFolder dfrFree, root:
	Duplicate/O outputMat, $(nameOfWave(inputMat)+ "_Ecorr")
	SumX24a($(nameOfWave(inputMat)+ "_Ecorr"))
endIf
End

Function getScofieldCrossSections()
NewPath/C/Q/Z/O IgorUserProc, (SpecialDirPath("Igor Pro User Files",0,0,0)+":Igor Pro 6 User Files:User Procedures")
Variable fileRefNum, p1, numPoints = 1, i = 0
String stringHolder, element, column_title_list, values_list,variables_list
DFREF savDF = GetDataFolderDFR(), csDF
NewDataFolder/O/S root:Scofield
csDF = GetDataFolderDFR()
Open/R/P=IgorUserProc/Z fileRefNum as "CrossSec_PE_Scofield.dat"
do
FReadLine fileRefNum, stringHolder
FStatus fileRefNum
if(char2num(stringHolder[0])==35)
	if(char2num(stringHolder[1])==83)
		SetDataFolder csDF
		sscanf stringHolder, "#S %*d %s", element
		NewDataFolder/O/S $(element)
		numPoints = 1
		i=0
	elseif(char2num(stringHolder[1])==76)
		column_title_list = ""
		For(p1 = 0; p1 < itemsinlist(stringHolder[2,inf], " "); p1 +=1)
			if(strlen(stringfromlist(p1, stringHolder[2,inf], " ")))
				if(strLen(stringfromlist(p1, stringHolder[2,inf], " "))!=0)
					column_title_list = AddListItem(replaceString("/",replaceString("\r",stringfromlist(p1, stringHolder[2,inf], " "), ""),""), column_title_list, ";", inf)
					Variable/G $(replaceString("/",replaceString("\r",stringfromlist(p1, stringHolder[2,inf], " "), ""),"") + "_Var")
					Make/O/N=(1) $(replaceString("/",replaceString("\r",stringfromlist(p1, StringHolder[2,inf], " "), ""),""))
				endif
			else
			endif
		endFor
	endif
else
	values_list = ""
	if (i >= numPoints)
		numPoints += 1
	endif
	For(p1 = 0; p1 < itemsinlist(stringHolder, " "); p1 += 1)
		if(strlen(stringfromlist(p1, stringHolder, " ")))
			values_list = AddListItem(replaceString("\r",stringfromlist(p1, stringHolder, " "), ""), values_list, ";", inf) 
		endif
	endFor
	For(p1 = 0; p1 < itemsInList(values_list); p1 += 1)
		NVAR holderVar = $(stringFromList(p1,column_title_list) + "_Var")
		Wave holderWave = $(stringFromList(p1,column_title_list))
		holderVar = str2num(stringFromList(p1,values_list))
		Redimension/N=(numPoints) holderWave
		holderWave[i] = holderVar
	endFor
	i += 1
endif
while(V_filePos < V_logEOF)
Close fileRefNum
SetDataFolder savDF
end

Menu "Spectrometry"
	"Load Scofield Cross Section", getScofieldCrossSections()
end

Function scaleTracesOnGraph(targetGraph)
String targetGraph
If(!strlen(targetGraph))
	targetGraph = WinName(0,1)
endIf
String listOfTraces = TraceNameList(targetGraph,";",1), listOfDupWaves = ""
Variable numberOfTraces = itemsInlist(listOfTraces), p1, waveVariance = inf, avgValue
Make/O/N=(numberOfTraces) waveVarianceWave
Make/O/N=(numberOfTraces)/WAVE waveReferenceWave

For(p1 = 0; p1 < numberOfTraces; p1 += 1)
	Wave tempWave = TraceNameToWaveRef(targetGraph,stringFromList(p1, listOfTraces))
	DFREF waveDFR =  GetWavesDataFolderDFR(tempWave)
	Duplicate/O tempWave, waveDFR:$(NameOfWave(tempWave) + "_N")
	waveReferenceWave[p1] = waveDFR:$(NameOfWave(tempWave) + "_N")
	Variable tempWaveVariance = Variance(tempWave,pnt2x(tempWave,numPnts(tempWave)),pnt2x(tempWave,numPnts(tempWave)-20))
	waveVarianceWave[p1] = tempWaveVariance
	If(abs(tempWaveVariance) < abs(waveVariance))
		waveVariance = tempWaveVariance
	endIf
endFor

Preferences 1;
For(p1 = 0; p1 < numberOfTraces; p1 += 1)
	Wave tempWave = waveReferenceWave[p1]
	tempWave *= sqrt(waveVariance/waveVarianceWave[p1])
	avgValue = Sum(tempWave, pnt2x(tempWave,numPnts(tempWave)),pnt2x(tempWave,numPnts(tempWave)-20))/20
	tempWave -=avgValue
	String traceName = stringFromList(p1, listOfTraces), xWave = XWaveName(targetGraph,traceName)
	if(strlen(xWave))
		xWave = GetWavesDataFolder(XWaveRefFromTrace(targetGraph,traceName),2)
	else
		xWave = ""
	endIf
	newGraphSameStyle(tempWave,xWave,targetGraph,traceName, typeOfGraph = 1)
endFor
Preferences 0;
KillWaves waveVarianceWave,waveReferenceWave
end